Reference documentation for deal.II version GIT 2c37b3e359 2022-05-26 00:10:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
grid_out.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
19 #include <deal.II/base/point.h>
22 
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)
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)
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)
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)
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)
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 
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.fail() == false, 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.fail() == false, ExcIO());
1013 }
1014 
1015 
1016 
1017 template <int dim, int spacedim>
1018 void
1020  std::ostream & out) const
1021 {
1022  AssertThrow(out.fail() == false, 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() +
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.fail() == false, ExcIO());
1114 }
1115 
1116 
1117 template <int dim, int spacedim>
1118 void
1120  std::ostream & out) const
1121 {
1122  AssertThrow(out.fail() == false, 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() +
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.fail() == false, 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.fail() == false, 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  // Now for each individual pair of dim and entry, add a physical group, if
1599  // necessary
1600  for (const auto &it : dim_entity_tag)
1601  {
1602  const auto &d = it.first;
1603  const auto &entity_tag = it.second;
1604  const auto &boundary_id = all_pairs[entity_tag - 1].first;
1605  const auto &manifold_id = all_pairs[entity_tag - 1].second;
1606 
1607  std::string physical_name;
1608  if (d == dim && boundary_id != 0)
1609  physical_name += "MaterialID:" + Utilities::int_to_string(
1610  static_cast<int>(boundary_id));
1611  else if (d < dim && boundary_id != 0)
1612  physical_name +=
1613  "BoundaryID:" +
1615  "-1" :
1616  Utilities::int_to_string(static_cast<int>(boundary_id)));
1617 
1618  std::string sep = physical_name != "" ? ", " : "";
1620  physical_name +=
1621  sep + "ManifoldID:" +
1622  Utilities::int_to_string(static_cast<int>(manifold_id));
1623  const auto physical_tag =
1624  gmsh::model::addPhysicalGroup(d, {entity_tag}, -1);
1625  if (physical_name != "")
1626  gmsh::model::setPhysicalName(d, physical_tag, physical_name);
1627  }
1628 
1629 
1630  gmsh::write(filename);
1631  gmsh::clear();
1632  gmsh::finalize();
1633 }
1634 #endif
1635 
1636 
1637 
1638 namespace
1639 {
1650  Point<2>
1651  svg_project_point(const Point<3> & point,
1652  const Point<3> & camera_position,
1653  const Tensor<1, 3> &camera_direction,
1654  const Tensor<1, 3> &camera_horizontal,
1655  const float camera_focus)
1656  {
1657  const Tensor<1, 3> camera_vertical =
1658  cross_product_3d(camera_horizontal, camera_direction);
1659 
1660  const float phi =
1661  camera_focus / ((point - camera_position) * camera_direction);
1662 
1663  const Point<3> projection =
1664  camera_position + phi * (point - camera_position);
1665 
1666  return {(projection - camera_position - camera_focus * camera_direction) *
1667  camera_horizontal,
1668  (projection - camera_position - camera_focus * camera_direction) *
1669  camera_vertical};
1670  }
1671 } // namespace
1672 
1673 
1674 
1675 template <int dim, int spacedim>
1676 void
1678  std::ostream & /*out*/) const
1679 {
1680  Assert(false,
1681  ExcMessage("Mesh output in SVG format is not implemented for anything "
1682  "other than two-dimensional meshes in two-dimensional "
1683  "space. That's because three-dimensional meshes are best "
1684  "viewed in programs that allow changing the viewpoint, "
1685  "but SVG format does not allow this: It is an inherently "
1686  "2d format, and for three-dimensional meshes would "
1687  "require choosing one, fixed viewpoint."
1688  "\n\n"
1689  "You probably want to output your mesh in a format such "
1690  "as VTK, VTU, or gnuplot."));
1691 }
1692 
1693 
1694 void
1695 GridOut::write_svg(const Triangulation<2, 2> &tria, std::ostream &out) const
1696 {
1697  unsigned int n = 0;
1698 
1699  unsigned int min_level, max_level;
1700 
1701  // Svg files require an underlying drawing grid. The size of this
1702  // grid is provided in the parameters height and width. Each of them
1703  // may be zero, such that it is computed from the other. Obviously,
1704  // both of them zero does not produce reasonable output.
1705  unsigned int height = svg_flags.height;
1706  unsigned int width = svg_flags.width;
1707  Assert(height != 0 || width != 0,
1708  ExcMessage("You have to set at least one of width and height"));
1709 
1710  unsigned int margin_in_percent = 0;
1712  margin_in_percent = 8;
1713 
1714  // initial font size for cell labels
1715  unsigned int cell_label_font_size;
1716 
1717  // get date and time
1718  // time_t time_stamp;
1719  // tm *now;
1720  // time_stamp = time(0);
1721  // now = localtime(&time_stamp);
1722 
1723  float camera_focus;
1724 
1725  Point<3> point;
1726  Point<2> projection_decomposition;
1727 
1728  float x_max_perspective, x_min_perspective;
1729  float y_max_perspective, y_min_perspective;
1730 
1731  float x_dimension_perspective, y_dimension_perspective;
1732 
1733 
1734  // auxiliary variables for the bounding box and the range of cell levels
1735  double x_min = tria.begin()->vertex(0)[0];
1736  double x_max = x_min;
1737  double y_min = tria.begin()->vertex(0)[1];
1738  double y_max = y_min;
1739 
1740  double x_dimension, y_dimension;
1741 
1742  min_level = max_level = tria.begin()->level();
1743 
1744  // auxiliary set for the materials being used
1745  std::set<unsigned int> materials;
1746 
1747  // auxiliary set for the levels being used
1748  std::set<unsigned int> levels;
1749 
1750  // auxiliary set for the subdomains being used
1751  std::set<unsigned int> subdomains;
1752 
1753  // auxiliary set for the level subdomains being used
1754  std::set<int> level_subdomains;
1755 
1756  // We use an active cell iterator to determine the
1757  // bounding box of the given triangulation and check
1758  // the cells for material id, level number, subdomain id
1759  // (, and level subdomain id).
1760  for (const auto &cell : tria.cell_iterators())
1761  {
1762  for (unsigned int vertex_index = 0; vertex_index < cell->n_vertices();
1763  ++vertex_index)
1764  {
1765  if (cell->vertex(vertex_index)[0] < x_min)
1766  x_min = cell->vertex(vertex_index)[0];
1767  if (cell->vertex(vertex_index)[0] > x_max)
1768  x_max = cell->vertex(vertex_index)[0];
1769 
1770  if (cell->vertex(vertex_index)[1] < y_min)
1771  y_min = cell->vertex(vertex_index)[1];
1772  if (cell->vertex(vertex_index)[1] > y_max)
1773  y_max = cell->vertex(vertex_index)[1];
1774  }
1775 
1776  if (static_cast<unsigned int>(cell->level()) < min_level)
1777  min_level = cell->level();
1778  if (static_cast<unsigned int>(cell->level()) > max_level)
1779  max_level = cell->level();
1780 
1781  materials.insert(cell->material_id());
1782  levels.insert(cell->level());
1783  if (cell->is_active())
1784  subdomains.insert(cell->subdomain_id() + 2);
1785  level_subdomains.insert(cell->level_subdomain_id() + 2);
1786  }
1787 
1788  x_dimension = x_max - x_min;
1789  y_dimension = y_max - y_min;
1790 
1791  // count the materials being used
1792  const unsigned int n_materials = materials.size();
1793 
1794  // count the levels being used
1795  const unsigned int n_levels = levels.size();
1796 
1797  // count the subdomains being used
1798  const unsigned int n_subdomains = subdomains.size();
1799 
1800  // count the level subdomains being used
1801  const unsigned int n_level_subdomains = level_subdomains.size();
1802 
1803  switch (svg_flags.coloring)
1804  {
1806  n = n_materials;
1807  break;
1809  n = n_levels;
1810  break;
1812  n = n_subdomains;
1813  break;
1815  n = n_level_subdomains;
1816  break;
1817  default:
1818  break;
1819  }
1820 
1821  // set the camera position to top view, targeting at the origin
1822  // vectors and variables for the perspective view
1823  Point<3> camera_position;
1824  camera_position[0] = 0;
1825  camera_position[1] = 0;
1826  camera_position[2] = 2. * std::max(x_dimension, y_dimension);
1827 
1828  Tensor<1, 3> camera_direction;
1829  camera_direction[0] = 0;
1830  camera_direction[1] = 0;
1831  camera_direction[2] = -1;
1832 
1833  Tensor<1, 3> camera_horizontal;
1834  camera_horizontal[0] = 1;
1835  camera_horizontal[1] = 0;
1836  camera_horizontal[2] = 0;
1837 
1838  camera_focus = .5 * std::max(x_dimension, y_dimension);
1839 
1840  Point<3> camera_position_temp;
1841  Point<3> camera_direction_temp;
1842  Point<3> camera_horizontal_temp;
1843 
1844  const double angle_factor = 3.14159265 / 180.;
1845 
1846  // (I) rotate the camera to the chosen polar angle
1847  camera_position_temp[1] =
1848  std::cos(angle_factor * svg_flags.polar_angle) * camera_position[1] -
1849  std::sin(angle_factor * svg_flags.polar_angle) * camera_position[2];
1850  camera_position_temp[2] =
1851  std::sin(angle_factor * svg_flags.polar_angle) * camera_position[1] +
1852  std::cos(angle_factor * svg_flags.polar_angle) * camera_position[2];
1853 
1854  camera_direction_temp[1] =
1855  std::cos(angle_factor * svg_flags.polar_angle) * camera_direction[1] -
1856  std::sin(angle_factor * svg_flags.polar_angle) * camera_direction[2];
1857  camera_direction_temp[2] =
1858  std::sin(angle_factor * svg_flags.polar_angle) * camera_direction[1] +
1859  std::cos(angle_factor * svg_flags.polar_angle) * camera_direction[2];
1860 
1861  camera_horizontal_temp[1] =
1862  std::cos(angle_factor * svg_flags.polar_angle) * camera_horizontal[1] -
1863  std::sin(angle_factor * svg_flags.polar_angle) * camera_horizontal[2];
1864  camera_horizontal_temp[2] =
1865  std::sin(angle_factor * svg_flags.polar_angle) * camera_horizontal[1] +
1866  std::cos(angle_factor * svg_flags.polar_angle) * camera_horizontal[2];
1867 
1868  camera_position[1] = camera_position_temp[1];
1869  camera_position[2] = camera_position_temp[2];
1870 
1871  camera_direction[1] = camera_direction_temp[1];
1872  camera_direction[2] = camera_direction_temp[2];
1873 
1874  camera_horizontal[1] = camera_horizontal_temp[1];
1875  camera_horizontal[2] = camera_horizontal_temp[2];
1876 
1877  // (II) rotate the camera to the chosen azimuth angle
1878  camera_position_temp[0] =
1879  std::cos(angle_factor * svg_flags.azimuth_angle) * camera_position[0] -
1880  std::sin(angle_factor * svg_flags.azimuth_angle) * camera_position[1];
1881  camera_position_temp[1] =
1882  std::sin(angle_factor * svg_flags.azimuth_angle) * camera_position[0] +
1883  std::cos(angle_factor * svg_flags.azimuth_angle) * camera_position[1];
1884 
1885  camera_direction_temp[0] =
1886  std::cos(angle_factor * svg_flags.azimuth_angle) * camera_direction[0] -
1887  std::sin(angle_factor * svg_flags.azimuth_angle) * camera_direction[1];
1888  camera_direction_temp[1] =
1889  std::sin(angle_factor * svg_flags.azimuth_angle) * camera_direction[0] +
1890  std::cos(angle_factor * svg_flags.azimuth_angle) * camera_direction[1];
1891 
1892  camera_horizontal_temp[0] =
1893  std::cos(angle_factor * svg_flags.azimuth_angle) * camera_horizontal[0] -
1894  std::sin(angle_factor * svg_flags.azimuth_angle) * camera_horizontal[1];
1895  camera_horizontal_temp[1] =
1896  std::sin(angle_factor * svg_flags.azimuth_angle) * camera_horizontal[0] +
1897  std::cos(angle_factor * svg_flags.azimuth_angle) * camera_horizontal[1];
1898 
1899  camera_position[0] = camera_position_temp[0];
1900  camera_position[1] = camera_position_temp[1];
1901 
1902  camera_direction[0] = camera_direction_temp[0];
1903  camera_direction[1] = camera_direction_temp[1];
1904 
1905  camera_horizontal[0] = camera_horizontal_temp[0];
1906  camera_horizontal[1] = camera_horizontal_temp[1];
1907 
1908  // translate the camera to the given triangulation
1909  camera_position[0] = x_min + .5 * x_dimension;
1910  camera_position[1] = y_min + .5 * y_dimension;
1911 
1912  camera_position[0] += 2. * std::max(x_dimension, y_dimension) *
1913  std::sin(angle_factor * svg_flags.polar_angle) *
1914  std::sin(angle_factor * svg_flags.azimuth_angle);
1915  camera_position[1] -= 2. * std::max(x_dimension, y_dimension) *
1916  std::sin(angle_factor * svg_flags.polar_angle) *
1917  std::cos(angle_factor * svg_flags.azimuth_angle);
1918 
1919 
1920  // determine the bounding box of the given triangulation on the projection
1921  // plane of the camera viewing system
1922  point[0] = tria.begin()->vertex(0)[0];
1923  point[1] = tria.begin()->vertex(0)[1];
1924  point[2] = 0;
1925 
1926  float min_level_min_vertex_distance = 0;
1927 
1929  {
1931  (static_cast<float>(tria.begin()->level()) /
1932  static_cast<float>(n_levels)) *
1933  std::max(x_dimension, y_dimension);
1934  }
1935 
1936  projection_decomposition = svg_project_point(
1937  point, camera_position, camera_direction, camera_horizontal, camera_focus);
1938 
1939  x_max_perspective = projection_decomposition[0];
1940  x_min_perspective = projection_decomposition[0];
1941 
1942  y_max_perspective = projection_decomposition[1];
1943  y_min_perspective = projection_decomposition[1];
1944 
1945  for (const auto &cell : tria.cell_iterators())
1946  {
1947  point[0] = cell->vertex(0)[0];
1948  point[1] = cell->vertex(0)[1];
1949  point[2] = 0;
1950 
1952  {
1953  point[2] =
1955  (static_cast<float>(cell->level()) / static_cast<float>(n_levels)) *
1956  std::max(x_dimension, y_dimension);
1957  }
1958 
1959  projection_decomposition = svg_project_point(point,
1960  camera_position,
1961  camera_direction,
1962  camera_horizontal,
1963  camera_focus);
1964 
1965  if (x_max_perspective < projection_decomposition[0])
1966  x_max_perspective = projection_decomposition[0];
1967  if (x_min_perspective > projection_decomposition[0])
1968  x_min_perspective = projection_decomposition[0];
1969 
1970  if (y_max_perspective < projection_decomposition[1])
1971  y_max_perspective = projection_decomposition[1];
1972  if (y_min_perspective > projection_decomposition[1])
1973  y_min_perspective = projection_decomposition[1];
1974 
1975  point[0] = cell->vertex(1)[0];
1976  point[1] = cell->vertex(1)[1];
1977 
1978  projection_decomposition = svg_project_point(point,
1979  camera_position,
1980  camera_direction,
1981  camera_horizontal,
1982  camera_focus);
1983 
1984  if (x_max_perspective < projection_decomposition[0])
1985  x_max_perspective = projection_decomposition[0];
1986  if (x_min_perspective > projection_decomposition[0])
1987  x_min_perspective = projection_decomposition[0];
1988 
1989  if (y_max_perspective < projection_decomposition[1])
1990  y_max_perspective = projection_decomposition[1];
1991  if (y_min_perspective > projection_decomposition[1])
1992  y_min_perspective = projection_decomposition[1];
1993 
1994  point[0] = cell->vertex(2)[0];
1995  point[1] = cell->vertex(2)[1];
1996 
1997  projection_decomposition = svg_project_point(point,
1998  camera_position,
1999  camera_direction,
2000  camera_horizontal,
2001  camera_focus);
2002 
2003  if (x_max_perspective < projection_decomposition[0])
2004  x_max_perspective = projection_decomposition[0];
2005  if (x_min_perspective > projection_decomposition[0])
2006  x_min_perspective = projection_decomposition[0];
2007 
2008  if (y_max_perspective < projection_decomposition[1])
2009  y_max_perspective = projection_decomposition[1];
2010  if (y_min_perspective > projection_decomposition[1])
2011  y_min_perspective = projection_decomposition[1];
2012 
2013  if (cell->n_vertices() == 4) // in case of quadrilateral
2014  {
2015  point[0] = cell->vertex(3)[0];
2016  point[1] = cell->vertex(3)[1];
2017 
2018  projection_decomposition = svg_project_point(point,
2019  camera_position,
2020  camera_direction,
2021  camera_horizontal,
2022  camera_focus);
2023 
2024  if (x_max_perspective < projection_decomposition[0])
2025  x_max_perspective = projection_decomposition[0];
2026  if (x_min_perspective > projection_decomposition[0])
2027  x_min_perspective = projection_decomposition[0];
2028 
2029  if (y_max_perspective < projection_decomposition[1])
2030  y_max_perspective = projection_decomposition[1];
2031  if (y_min_perspective > projection_decomposition[1])
2032  y_min_perspective = projection_decomposition[1];
2033  }
2034 
2035  if (static_cast<unsigned int>(cell->level()) == min_level)
2036  min_level_min_vertex_distance = cell->minimum_vertex_distance();
2037  }
2038 
2039  x_dimension_perspective = x_max_perspective - x_min_perspective;
2040  y_dimension_perspective = y_max_perspective - y_min_perspective;
2041 
2042  // create the svg file with an internal style sheet
2043  if (width == 0)
2044  width = static_cast<unsigned int>(
2045  .5 + height * (x_dimension_perspective / y_dimension_perspective));
2046  else if (height == 0)
2047  height = static_cast<unsigned int>(
2048  .5 + width * (y_dimension_perspective / x_dimension_perspective));
2049  unsigned int additional_width = 0;
2050  // font size for date, time, legend, and colorbar
2051  unsigned int font_size =
2052  static_cast<unsigned int>(.5 + (height / 100.) * 1.75);
2053  cell_label_font_size = static_cast<unsigned int>(
2054  .5 + (height * .15 * svg_flags.cell_font_scaling *
2055  min_level_min_vertex_distance / std::min(x_dimension, y_dimension)));
2056 
2057  if (svg_flags.draw_legend &&
2061  {
2062  additional_width = static_cast<unsigned int>(
2063  .5 + height * .4); // additional width for legend
2064  }
2065  else if (svg_flags.draw_colorbar && (svg_flags.coloring != 0u))
2066  {
2067  additional_width = static_cast<unsigned int>(
2068  .5 + height * .175); // additional width for colorbar
2069  }
2070 
2071  // out << "<!-- deal.ii GridOut " << now->tm_mday << '/' << now->tm_mon + 1 <<
2072  // '/' << now->tm_year + 1900
2073  // << ' ' << now->tm_hour << ':';
2074  //
2075  // if (now->tm_min < 10) out << '0';
2076  //
2077  // out << now->tm_min << " -->" << '\n';
2078 
2079  // basic svg header
2080  out << "<svg width=\"" << width + additional_width << "\" height=\"" << height
2081  << "\" xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\">" << '\n'
2082  << '\n';
2083 
2084 
2086  {
2087  out
2088  << " <linearGradient id=\"background_gradient\" gradientUnits=\"userSpaceOnUse\" x1=\"0\" y1=\"0\" x2=\"0\" y2=\""
2089  << height << "\">" << '\n'
2090  << " <stop offset=\"0\" style=\"stop-color:white\"/>" << '\n'
2091  << " <stop offset=\"1\" style=\"stop-color:lightsteelblue\"/>" << '\n'
2092  << " </linearGradient>" << '\n';
2093  }
2094 
2095  out << '\n';
2096 
2097  // header for the internal style sheet
2098  out << "<!-- internal style sheet -->" << '\n'
2099  << "<style type=\"text/css\"><![CDATA[" << '\n';
2100 
2101  // set the background of the output graphic
2103  out << " rect.background{fill:url(#background_gradient)}" << '\n';
2105  out << " rect.background{fill:white}" << '\n';
2106  else
2107  out << " rect.background{fill:none}" << '\n';
2108 
2109  // basic svg graphic element styles
2110  out << " rect{fill:none; stroke:rgb(25,25,25); stroke-width:"
2111  << svg_flags.line_thickness << '}' << '\n'
2112  << " text{font-family:Helvetica; text-anchor:middle; fill:rgb(25,25,25)}"
2113  << '\n'
2114  << " line{stroke:rgb(25,25,25); stroke-width:"
2115  << svg_flags.boundary_line_thickness << '}' << '\n'
2116  << " path{fill:none; stroke:rgb(25,25,25); stroke-width:"
2117  << svg_flags.line_thickness << '}' << '\n'
2118  << " circle{fill:white; stroke:black; stroke-width:2}" << '\n'
2119  << '\n';
2120 
2121  // polygon styles with respect to the chosen cell coloring
2122  if (svg_flags.coloring != 0u)
2123  {
2124  unsigned int labeling_index = 0;
2125  auto materials_it = materials.begin();
2126  auto levels_it = levels.begin();
2127  auto subdomains_it = subdomains.begin();
2128  auto level_subdomains_it = level_subdomains.begin();
2129 
2130  for (unsigned int index = 0; index < n; ++index)
2131  {
2132  double h;
2133 
2134  if (n != 1)
2135  h = .6 - (index / (n - 1.)) * .6;
2136  else
2137  h = .6;
2138 
2139  unsigned int r = 0;
2140  unsigned int g = 0;
2141  unsigned int b = 0;
2142 
2143  unsigned int i = static_cast<unsigned int>(h * 6);
2144 
2145  double f = h * 6 - i;
2146  double q = 1 - f;
2147  double t = f;
2148 
2149  switch (i % 6)
2150  {
2151  case 0:
2152  r = 255, g = static_cast<unsigned int>(.5 + 255 * t);
2153  break;
2154  case 1:
2155  r = static_cast<unsigned int>(.5 + 255 * q), g = 255;
2156  break;
2157  case 2:
2158  g = 255, b = static_cast<unsigned int>(.5 + 255 * t);
2159  break;
2160  case 3:
2161  g = static_cast<unsigned int>(.5 + 255 * q), b = 255;
2162  break;
2163  case 4:
2164  r = static_cast<unsigned int>(.5 + 255 * t), b = 255;
2165  break;
2166  case 5:
2167  r = 255, b = static_cast<unsigned int>(.5 + 255 * q);
2168  break;
2169  default:
2170  break;
2171  }
2172 
2173  switch (svg_flags.coloring)
2174  {
2176  labeling_index = *materials_it++;
2177  break;
2179  labeling_index = *levels_it++;
2180  break;
2182  labeling_index = *subdomains_it++;
2183  break;
2185  labeling_index = *level_subdomains_it++;
2186  break;
2187  default:
2188  break;
2189  }
2190 
2191  out << " path.p" << labeling_index << "{fill:rgb(" << r << ',' << g
2192  << ',' << b << "); "
2193  << "stroke:rgb(25,25,25); stroke-width:"
2194  << svg_flags.line_thickness << '}' << '\n';
2195 
2196  out << " path.ps" << labeling_index << "{fill:rgb("
2197  << static_cast<unsigned int>(.5 + .75 * r) << ','
2198  << static_cast<unsigned int>(.5 + .75 * g) << ','
2199  << static_cast<unsigned int>(.5 + .75 * b) << "); "
2200  << "stroke:rgb(20,20,20); stroke-width:"
2201  << svg_flags.line_thickness << '}' << '\n';
2202 
2203  out << " rect.r" << labeling_index << "{fill:rgb(" << r << ',' << g
2204  << ',' << b << "); "
2205  << "stroke:rgb(25,25,25); stroke-width:"
2206  << svg_flags.line_thickness << '}' << '\n';
2207 
2208  labeling_index++;
2209  }
2210  }
2211 
2212  out << "]]></style>" << '\n' << '\n';
2213 
2214  // background rectangle
2215  out << " <rect class=\"background\" width=\"" << width << "\" height=\""
2216  << height << "\"/>" << '\n';
2217 
2219  {
2220  unsigned int x_offset = 0;
2221 
2222  if (svg_flags.margin)
2223  x_offset = static_cast<unsigned int>(.5 + (height / 100.) *
2224  (margin_in_percent / 2.));
2225  else
2226  x_offset = static_cast<unsigned int>(.5 + height * .025);
2227 
2228  out
2229  << " <text x=\"" << x_offset << "\" y=\""
2230  << static_cast<unsigned int>(.5 + height * .0525) << '\"'
2231  << " style=\"font-weight:100; fill:lightsteelblue; text-anchor:start; font-family:Courier; font-size:"
2232  << static_cast<unsigned int>(.5 + height * .045) << "px\">"
2233  << "deal.II"
2234  << "</text>" << '\n';
2235 
2236  // out << " <text x=\"" << x_offset + static_cast<unsigned int>(.5 +
2237  // height * .045 * 4.75) << "\" y=\"" << static_cast<unsigned int>(.5 +
2238  // height * .0525) << '\"'
2239  // << " style=\"fill:lightsteelblue; text-anchor:start; font-size:" <<
2240  // font_size << "\">"
2241  // << now->tm_mday << '/' << now->tm_mon + 1 << '/' << now->tm_year +
2242  // 1900
2243  // << " - " << now->tm_hour << ':';
2244  //
2245  // if(now->tm_min < 10) out << '0';
2246  //
2247  // out << now->tm_min
2248  // << "</text>"<< '\n' << '\n';
2249  }
2250 
2251  // draw the cells, starting out from the minimal level (in order to guaranty a
2252  // correct perspective view)
2253  out << " <!-- cells -->" << '\n';
2254 
2255  for (unsigned int level_index = min_level; level_index <= max_level;
2256  level_index++)
2257  {
2258  for (const auto &cell : tria.cell_iterators_on_level(level_index))
2259  {
2260  if (!svg_flags.convert_level_number_to_height && !cell->is_active())
2261  continue;
2262 
2263  // draw the current cell
2264  out << " <path";
2265 
2266  if (svg_flags.coloring != 0u)
2267  {
2268  out << " class=\"p";
2269 
2270  if (!cell->is_active() &&
2272  out << 's';
2273 
2274  switch (svg_flags.coloring)
2275  {
2277  out << cell->material_id();
2278  break;
2280  out << static_cast<unsigned int>(cell->level());
2281  break;
2283  if (cell->is_active())
2284  out << cell->subdomain_id() + 2;
2285  else
2286  out << 'X';
2287  break;
2289  out << cell->level_subdomain_id() + 2;
2290  break;
2291  default:
2292  break;
2293  }
2294 
2295  out << '\"';
2296  }
2297 
2298  out << " d=\"M ";
2299 
2300  point[0] = cell->vertex(0)[0];
2301  point[1] = cell->vertex(0)[1];
2302  point[2] = 0;
2303 
2305  {
2307  (static_cast<float>(cell->level()) /
2308  static_cast<float>(n_levels)) *
2309  std::max(x_dimension, y_dimension);
2310  }
2311 
2312  projection_decomposition = svg_project_point(point,
2313  camera_position,
2314  camera_direction,
2315  camera_horizontal,
2316  camera_focus);
2317 
2318  out << static_cast<unsigned int>(
2319  .5 +
2320  ((projection_decomposition[0] - x_min_perspective) /
2321  x_dimension_perspective) *
2322  (width - (width / 100.) * 2. * margin_in_percent) +
2323  ((width / 100.) * margin_in_percent))
2324  << ' '
2325  << static_cast<unsigned int>(
2326  .5 + height - (height / 100.) * margin_in_percent -
2327  ((projection_decomposition[1] - y_min_perspective) /
2328  y_dimension_perspective) *
2329  (height - (height / 100.) * 2. * margin_in_percent));
2330 
2331  out << " L ";
2332 
2333  point[0] = cell->vertex(1)[0];
2334  point[1] = cell->vertex(1)[1];
2335 
2336  projection_decomposition = svg_project_point(point,
2337  camera_position,
2338  camera_direction,
2339  camera_horizontal,
2340  camera_focus);
2341 
2342  out << static_cast<unsigned int>(
2343  .5 +
2344  ((projection_decomposition[0] - x_min_perspective) /
2345  x_dimension_perspective) *
2346  (width - (width / 100.) * 2. * margin_in_percent) +
2347  ((width / 100.) * margin_in_percent))
2348  << ' '
2349  << static_cast<unsigned int>(
2350  .5 + height - (height / 100.) * margin_in_percent -
2351  ((projection_decomposition[1] - y_min_perspective) /
2352  y_dimension_perspective) *
2353  (height - (height / 100.) * 2. * margin_in_percent));
2354 
2355  out << " L ";
2356 
2357  if (cell->n_vertices() == 4) // in case of quadrilateral
2358  {
2359  point[0] = cell->vertex(3)[0];
2360  point[1] = cell->vertex(3)[1];
2361 
2362  projection_decomposition = svg_project_point(point,
2363  camera_position,
2364  camera_direction,
2365  camera_horizontal,
2366  camera_focus);
2367 
2368  out << static_cast<unsigned int>(
2369  .5 +
2370  ((projection_decomposition[0] - x_min_perspective) /
2371  x_dimension_perspective) *
2372  (width - (width / 100.) * 2. * margin_in_percent) +
2373  ((width / 100.) * margin_in_percent))
2374  << ' '
2375  << static_cast<unsigned int>(
2376  .5 + height - (height / 100.) * margin_in_percent -
2377  ((projection_decomposition[1] - y_min_perspective) /
2378  y_dimension_perspective) *
2379  (height - (height / 100.) * 2. * margin_in_percent));
2380 
2381  out << " L ";
2382  }
2383 
2384  point[0] = cell->vertex(2)[0];
2385  point[1] = cell->vertex(2)[1];
2386 
2387  projection_decomposition = svg_project_point(point,
2388  camera_position,
2389  camera_direction,
2390  camera_horizontal,
2391  camera_focus);
2392 
2393  out << static_cast<unsigned int>(
2394  .5 +
2395  ((projection_decomposition[0] - x_min_perspective) /
2396  x_dimension_perspective) *
2397  (width - (width / 100.) * 2. * margin_in_percent) +
2398  ((width / 100.) * margin_in_percent))
2399  << ' '
2400  << static_cast<unsigned int>(
2401  .5 + height - (height / 100.) * margin_in_percent -
2402  ((projection_decomposition[1] - y_min_perspective) /
2403  y_dimension_perspective) *
2404  (height - (height / 100.) * 2. * margin_in_percent));
2405 
2406  out << " L ";
2407 
2408  point[0] = cell->vertex(0)[0];
2409  point[1] = cell->vertex(0)[1];
2410 
2411  projection_decomposition = svg_project_point(point,
2412  camera_position,
2413  camera_direction,
2414  camera_horizontal,
2415  camera_focus);
2416 
2417  out << static_cast<unsigned int>(
2418  .5 +
2419  ((projection_decomposition[0] - x_min_perspective) /
2420  x_dimension_perspective) *
2421  (width - (width / 100.) * 2. * margin_in_percent) +
2422  ((width / 100.) * margin_in_percent))
2423  << ' '
2424  << static_cast<unsigned int>(
2425  .5 + height - (height / 100.) * margin_in_percent -
2426  ((projection_decomposition[1] - y_min_perspective) /
2427  y_dimension_perspective) *
2428  (height - (height / 100.) * 2. * margin_in_percent));
2429 
2430  out << "\"/>" << '\n';
2431 
2432  // label the current cell
2436  {
2437  point[0] = cell->center()[0];
2438  point[1] = cell->center()[1];
2439  point[2] = 0;
2440 
2442  {
2444  (static_cast<float>(cell->level()) /
2445  static_cast<float>(n_levels)) *
2446  std::max(x_dimension, y_dimension);
2447  }
2448 
2449  const double distance_to_camera =
2450  std::sqrt(std::pow(point[0] - camera_position[0], 2.) +
2451  std::pow(point[1] - camera_position[1], 2.) +
2452  std::pow(point[2] - camera_position[2], 2.));
2453  const double distance_factor =
2454  distance_to_camera / (2. * std::max(x_dimension, y_dimension));
2455 
2456  projection_decomposition = svg_project_point(point,
2457  camera_position,
2458  camera_direction,
2459  camera_horizontal,
2460  camera_focus);
2461 
2462  const unsigned int font_size_this_cell =
2463  static_cast<unsigned int>(
2464  .5 +
2465  cell_label_font_size *
2466  std::pow(.5, cell->level() - 4. + 3.5 * distance_factor));
2467 
2468  out << " <text"
2469  << " x=\""
2470  << static_cast<unsigned int>(
2471  .5 +
2472  ((projection_decomposition[0] - x_min_perspective) /
2473  x_dimension_perspective) *
2474  (width - (width / 100.) * 2. * margin_in_percent) +
2475  ((width / 100.) * margin_in_percent))
2476  << "\" y=\""
2477  << static_cast<unsigned int>(
2478  .5 + height - (height / 100.) * margin_in_percent -
2479  ((projection_decomposition[1] - y_min_perspective) /
2480  y_dimension_perspective) *
2481  (height - (height / 100.) * 2. * margin_in_percent) +
2482  0.5 * font_size_this_cell)
2483  << "\" style=\"font-size:" << font_size_this_cell << "px\">";
2484 
2486  {
2487  out << cell->level();
2488  }
2489 
2491  {
2493  out << '.';
2494  out << cell->index();
2495  }
2496 
2498  {
2501  out << ',';
2502  out
2503  << static_cast<std::make_signed<types::material_id>::type>(
2504  cell->material_id());
2505  }
2506 
2508  {
2511  out << ',';
2512  if (cell->is_active())
2513  out << static_cast<
2514  std::make_signed<types::subdomain_id>::type>(
2515  cell->subdomain_id());
2516  else
2517  out << 'X';
2518  }
2519 
2521  {
2526  out << ',';
2527  out
2528  << static_cast<std::make_signed<types::subdomain_id>::type>(
2529  cell->level_subdomain_id());
2530  }
2531 
2532  out << "</text>" << '\n';
2533  }
2534 
2535  // if the current cell lies at the boundary of the triangulation, draw
2536  // the additional boundary line
2538  {
2539  for (auto faceIndex : cell->face_indices())
2540  {
2541  if (cell->at_boundary(faceIndex))
2542  {
2543  point[0] = cell->face(faceIndex)->vertex(0)[0];
2544  point[1] = cell->face(faceIndex)->vertex(0)[1];
2545  point[2] = 0;
2546 
2548  {
2550  (static_cast<float>(cell->level()) /
2551  static_cast<float>(n_levels)) *
2552  std::max(x_dimension, y_dimension);
2553  }
2554 
2555  projection_decomposition =
2556  svg_project_point(point,
2557  camera_position,
2558  camera_direction,
2559  camera_horizontal,
2560  camera_focus);
2561 
2562  out << " <line x1=\""
2563  << static_cast<unsigned int>(
2564  .5 +
2565  ((projection_decomposition[0] -
2566  x_min_perspective) /
2567  x_dimension_perspective) *
2568  (width -
2569  (width / 100.) * 2. * margin_in_percent) +
2570  ((width / 100.) * margin_in_percent))
2571  << "\" y1=\""
2572  << static_cast<unsigned int>(
2573  .5 + height -
2574  (height / 100.) * margin_in_percent -
2575  ((projection_decomposition[1] -
2576  y_min_perspective) /
2577  y_dimension_perspective) *
2578  (height -
2579  (height / 100.) * 2. * margin_in_percent));
2580 
2581  point[0] = cell->face(faceIndex)->vertex(1)[0];
2582  point[1] = cell->face(faceIndex)->vertex(1)[1];
2583  point[2] = 0;
2584 
2586  {
2588  (static_cast<float>(cell->level()) /
2589  static_cast<float>(n_levels)) *
2590  std::max(x_dimension, y_dimension);
2591  }
2592 
2593  projection_decomposition =
2594  svg_project_point(point,
2595  camera_position,
2596  camera_direction,
2597  camera_horizontal,
2598  camera_focus);
2599 
2600  out << "\" x2=\""
2601  << static_cast<unsigned int>(
2602  .5 +
2603  ((projection_decomposition[0] -
2604  x_min_perspective) /
2605  x_dimension_perspective) *
2606  (width -
2607  (width / 100.) * 2. * margin_in_percent) +
2608  ((width / 100.) * margin_in_percent))
2609  << "\" y2=\""
2610  << static_cast<unsigned int>(
2611  .5 + height -
2612  (height / 100.) * margin_in_percent -
2613  ((projection_decomposition[1] -
2614  y_min_perspective) /
2615  y_dimension_perspective) *
2616  (height -
2617  (height / 100.) * 2. * margin_in_percent))
2618  << "\"/>" << '\n';
2619 
2620 
2622  {
2623  const double distance_to_camera = std::sqrt(
2624  std::pow(point[0] - camera_position[0], 2.) +
2625  std::pow(point[1] - camera_position[1], 2.) +
2626  std::pow(point[2] - camera_position[2], 2.));
2627  const double distance_factor =
2628  distance_to_camera /
2629  (2. * std::max(x_dimension, y_dimension));
2630 
2631  const unsigned int font_size_this_edge =
2632  static_cast<unsigned int>(
2633  .5 + .5 * cell_label_font_size *
2634  std::pow(.5,
2635  cell->level() - 4. +
2636  3.5 * distance_factor));
2637 
2638  point[0] = cell->face(faceIndex)->center()[0];
2639  point[1] = cell->face(faceIndex)->center()[1];
2640  point[2] = 0;
2641 
2643  {
2645  (static_cast<float>(cell->level()) /
2646  static_cast<float>(n_levels)) *
2647  std::max(x_dimension, y_dimension);
2648  }
2649 
2650  projection_decomposition =
2651  svg_project_point(point,
2652  camera_position,
2653  camera_direction,
2654  camera_horizontal,
2655  camera_focus);
2656 
2657  const unsigned int xc = static_cast<unsigned int>(
2658  .5 +
2659  ((projection_decomposition[0] - x_min_perspective) /
2660  x_dimension_perspective) *
2661  (width -
2662  (width / 100.) * 2. * margin_in_percent) +
2663  ((width / 100.) * margin_in_percent));
2664  const unsigned int yc = static_cast<unsigned int>(
2665  .5 + height - (height / 100.) * margin_in_percent -
2666  ((projection_decomposition[1] - y_min_perspective) /
2667  y_dimension_perspective) *
2668  (height -
2669  (height / 100.) * 2. * margin_in_percent));
2670 
2671  out << " <circle cx=\"" << xc << "\" cy=\"" << yc
2672  << "\" r=\"" << font_size_this_edge << "\" />"
2673  << '\n';
2674 
2675  out << " <text x=\"" << xc << "\" y=\"" << yc
2676  << "\" style=\"font-size:" << font_size_this_edge
2677  << "px\" dominant-baseline=\"middle\">"
2678  << static_cast<int>(
2679  cell->face(faceIndex)->boundary_id())
2680  << "</text>" << '\n';
2681  }
2682  }
2683  }
2684  }
2685  }
2686  }
2687 
2688 
2689 
2690  // draw the legend
2691  if (svg_flags.draw_legend)
2692  out << '\n' << " <!-- legend -->" << '\n';
2693 
2694  additional_width = 0;
2695  if (!svg_flags.margin)
2696  additional_width = static_cast<unsigned int>(.5 + (height / 100.) * 2.5);
2697 
2698  // explanation of the cell labeling
2699  if (svg_flags.draw_legend &&
2703  {
2704  unsigned int line_offset = 0;
2705  out << " <rect x=\"" << width + additional_width << "\" y=\""
2706  << static_cast<unsigned int>(.5 + (height / 100.) * margin_in_percent)
2707  << "\" width=\""
2708  << static_cast<unsigned int>(.5 + (height / 100.) *
2709  (40. - margin_in_percent))
2710  << "\" height=\"" << static_cast<unsigned int>(.5 + height * .215)
2711  << "\"/>" << '\n';
2712 
2713  out << " <text x=\""
2714  << width + additional_width +
2715  static_cast<unsigned int>(.5 + (height / 100.) * 1.25)
2716  << "\" y=\""
2717  << static_cast<unsigned int>(.5 +
2718  (height / 100.) * margin_in_percent +
2719  (++line_offset) * 1.5 * font_size)
2720  << "\" style=\"text-anchor:start; font-weight:bold; font-size:"
2721  << font_size << "px\">"
2722  << "cell label"
2723  << "</text>" << '\n';
2724 
2726  {
2727  out << " <text x=\""
2728  << width + additional_width +
2729  static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2730  << "\" y=\""
2731  << static_cast<unsigned int>(.5 +
2732  (height / 100.) * margin_in_percent +
2733  (++line_offset) * 1.5 * font_size)
2734  << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2735  << font_size << "px\">"
2736  << "cell_level";
2737 
2741  out << '.';
2742 
2743  out << "</text>" << '\n';
2744  }
2745 
2747  {
2748  out << " <text x=\""
2749  << width + additional_width +
2750  static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2751  << "\" y=\""
2752  << static_cast<unsigned int>(.5 +
2753  (height / 100.) * margin_in_percent +
2754  (++line_offset) * 1.5 * font_size)
2755  << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2756  << font_size << "px\">"
2757  << "cell_index";
2758 
2761  out << ',';
2762 
2763  out << "</text>" << '\n';
2764  }
2765 
2767  {
2768  out << " <text x=\""
2769  << width + additional_width +
2770  static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2771  << "\" y=\""
2772  << static_cast<unsigned int>(.5 +
2773  (height / 100.) * margin_in_percent +
2774  (++line_offset) * 1.5 * font_size)
2775  << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2776  << font_size << "px\">"
2777  << "material_id";
2778 
2781  out << ',';
2782 
2783  out << "</text>" << '\n';
2784  }
2785 
2787  {
2788  out << " <text x= \""
2789  << width + additional_width +
2790  static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2791  << "\" y=\""
2792  << static_cast<unsigned int>(.5 +
2793  (height / 100.) * margin_in_percent +
2794  (++line_offset) * 1.5 * font_size)
2795  << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2796  << font_size << "px\">"
2797  << "subdomain_id";
2798 
2800  out << ',';
2801 
2802  out << "</text>" << '\n';
2803  }
2804 
2806  {
2807  out << " <text x= \""
2808  << width + additional_width +
2809  static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2810  << "\" y=\""
2811  << static_cast<unsigned int>(.5 +
2812  (height / 100.) * margin_in_percent +
2813  (++line_offset) * 1.5 * font_size)
2814  << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2815  << font_size << "px\">"
2816  << "level_subdomain_id"
2817  << "</text>" << '\n';
2818  }
2819 
2821  {
2822  out << " <text x=\""
2823  << width + additional_width +
2824  static_cast<unsigned int>(.5 + (height / 100.) * 1.25)
2825  << "\" y=\""
2826  << static_cast<unsigned int>(.5 +
2827  (height / 100.) * margin_in_percent +
2828  (++line_offset) * 1.5 * font_size)
2829  << "\" style=\"text-anchor:start; font-weight:bold; font-size:"
2830  << font_size << "px\">"
2831  << "edge label"
2832  << "</text>" << '\n';
2833 
2834  out << " <text x= \""
2835  << width + additional_width +
2836  static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2837  << "\" y=\""
2838  << static_cast<unsigned int>(.5 +
2839  (height / 100.) * margin_in_percent +
2840  (++line_offset) * 1.5 * font_size)
2841  << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2842  << font_size << "px\">"
2843  << "boundary_id"
2844  << "</text>" << '\n';
2845  }
2846  }
2847 
2848  // show azimuth angle and polar angle as text below the explanation of the
2849  // cell labeling
2850  if (svg_flags.draw_legend)
2851  {
2852  out << " <text x=\"" << width + additional_width << "\" y=\""
2853  << static_cast<unsigned int>(
2854  .5 + (height / 100.) * margin_in_percent + 13.75 * font_size)
2855  << "\" style=\"text-anchor:start; font-size:" << font_size << "px\">"
2856  << "azimuth: " << svg_flags.azimuth_angle
2857  << "°, polar: " << svg_flags.polar_angle << "°</text>" << '\n';
2858  }
2859 
2860 
2861  // draw the colorbar
2862  if (svg_flags.draw_colorbar && (svg_flags.coloring != 0u))
2863  {
2864  out << '\n' << " <!-- colorbar -->" << '\n';
2865 
2866  out << " <text x=\"" << width + additional_width << "\" y=\""
2867  << static_cast<unsigned int>(
2868  .5 + (height / 100.) * (margin_in_percent + 29.) -
2869  (font_size / 1.25))
2870  << "\" style=\"text-anchor:start; font-weight:bold; font-size:"
2871  << font_size << "px\">";
2872 
2873  switch (svg_flags.coloring)
2874  {
2875  case 1:
2876  out << "material_id";
2877  break;
2878  case 2:
2879  out << "level_number";
2880  break;
2881  case 3:
2882  out << "subdomain_id";
2883  break;
2884  case 4:
2885  out << "level_subdomain_id";
2886  break;
2887  default:
2888  break;
2889  }
2890 
2891  out << "</text>" << '\n';
2892 
2893  unsigned int element_height = static_cast<unsigned int>(
2894  ((height / 100.) * (71. - 2. * margin_in_percent)) / n);
2895  unsigned int element_width =
2896  static_cast<unsigned int>(.5 + (height / 100.) * 2.5);
2897 
2898  int labeling_index = 0;
2899  auto materials_it = materials.begin();
2900  auto levels_it = levels.begin();
2901  auto subdomains_it = subdomains.begin();
2902  auto level_subdomains_it = level_subdomains.begin();
2903 
2904  for (unsigned int index = 0; index < n; ++index)
2905  {
2906  switch (svg_flags.coloring)
2907  {
2909  labeling_index = *materials_it++;
2910  break;
2912  labeling_index = *levels_it++;
2913  break;
2915  labeling_index = *subdomains_it++;
2916  break;
2918  labeling_index = *level_subdomains_it++;
2919  break;
2920  default:
2921  break;
2922  }
2923 
2924  out << " <rect class=\"r" << labeling_index << "\" x=\""
2925  << width + additional_width << "\" y=\""
2926  << static_cast<unsigned int>(.5 + (height / 100.) *
2927  (margin_in_percent + 29)) +
2928  (n - index - 1) * element_height
2929  << "\" width=\"" << element_width << "\" height=\""
2930  << element_height << "\"/>" << '\n';
2931 
2932  out << " <text x=\""
2933  << width + additional_width + 1.5 * element_width << "\" y=\""
2934  << static_cast<unsigned int>(.5 + (height / 100.) *
2935  (margin_in_percent + 29)) +
2936  (n - index - 1 + .5) * element_height +
2937  static_cast<unsigned int>(.5 + font_size * .35)
2938  << "\""
2939  << " style=\"text-anchor:start; font-size:"
2940  << static_cast<unsigned int>(.5 + font_size) << "px";
2941 
2942  if (index == 0 || index == n - 1)
2943  out << "; font-weight:bold";
2944 
2945  out << "\">" << labeling_index;
2946 
2947  if (index == n - 1)
2948  out << " max";
2949  if (index == 0)
2950  out << " min";
2951 
2952  out << "</text>" << '\n';
2953 
2954  labeling_index++;
2955  }
2956  }
2957 
2958 
2959  // finalize the svg file
2960  out << '\n' << "</svg>";
2961  out.flush();
2962 }
2963 
2964 
2965 
2966 template <>
2967 void
2968 GridOut::write_mathgl(const Triangulation<1> &, std::ostream &) const
2969 {
2970  // 1d specialization not done yet
2971  Assert(false, ExcNotImplemented());
2972 }
2973 
2974 
2975 
2976 template <int dim, int spacedim>
2977 void
2979  std::ostream & out) const
2980 {
2981  AssertThrow(out.fail() == false, ExcIO());
2982 
2983  // (i) write header
2984  {
2985  // block this to have local variables destroyed after use
2986  const std::time_t time1 = std::time(nullptr);
2987  const std::tm * time = std::localtime(&time1);
2988 
2989  out
2990  << "\n#"
2991  << "\n# This file was generated by the deal.II library."
2992  << "\n# Date = " << time->tm_year + 1900 << "/" << std::setfill('0')
2993  << std::setw(2) << time->tm_mon + 1 << "/" << std::setfill('0')
2994  << std::setw(2) << time->tm_mday << "\n# Time = " << std::setfill('0')
2995  << std::setw(2) << time->tm_hour << ":" << std::setfill('0')
2996  << std::setw(2) << time->tm_min << ":" << std::setfill('0')
2997  << std::setw(2) << time->tm_sec << "\n#"
2998  << "\n# For a description of the MathGL script format see the MathGL manual. "
2999  << "\n#"
3000  << "\n# Note: This file is understood by MathGL v2.1 and higher only, and can "
3001  << "\n# be quickly viewed in a graphical environment using \'mglview\'. "
3002  << "\n#" << '\n';
3003  }
3004 
3005  // define a helper to keep loops approximately dim-independent
3006  // since MathGL labels axes as x, y, z
3007  const std::string axes = "xyz";
3008 
3009  // (ii) write preamble and graphing tweaks
3010  out << "\n#"
3011  << "\n# Preamble."
3012  << "\n#" << '\n';
3013 
3015  out << "\nbox";
3016 
3017  // deal with dimension dependent preamble; eg. default sizes and
3018  // views for MathGL (cf. gnuplot).
3019  switch (dim)
3020  {
3021  case 2:
3022  out << "\nsetsize 800 800";
3023  out << "\nrotate 0 0";
3024  break;
3025  case 3:
3026  out << "\nsetsize 800 800";
3027  out << "\nrotate 60 40";
3028  break;
3029  default:
3030  Assert(false, ExcNotImplemented());
3031  }
3032  out << '\n';
3033 
3034  // (iii) write vertex ordering
3035  out << "\n#"
3036  << "\n# Vertex ordering."
3037  << "\n# list <vertex order> <vertex indices>"
3038  << "\n#" << '\n';
3039 
3040  // todo: This denotes the natural ordering of vertices, but it needs
3041  // to check this is really always true for a given grid (it's not
3042  // true in @ref step_1 "step-1" grid-2 for instance).
3043  switch (dim)
3044  {
3045  case 2:
3046  out << "\nlist f 0 1 2 3" << '\n';
3047  break;
3048  case 3:
3049  out
3050  << "\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"
3051  << '\n';
3052  break;
3053  default:
3054  Assert(false, ExcNotImplemented());
3055  }
3056 
3057  // (iv) write a list of vertices of cells
3058  out << "\n#"
3059  << "\n# List of vertices."
3060  << "\n# list <id> <vertices>"
3061  << "\n#" << '\n';
3062 
3063  // run over all active cells and write out a list of
3064  // xyz-coordinates that correspond to vertices
3065  // No global indices in deal.II, so we make one up here.
3066  for (const auto &cell : tria.active_cell_iterators())
3067  {
3068  for (unsigned int i = 0; i < dim; ++i)
3069  {
3070  // if (cell->direction_flag ()==true)
3071  // out << "\ntrue";
3072  // else
3073  // out << "\nfalse";
3074 
3075  out << "\nlist " << axes[i] << cell->active_cell_index() << " ";
3076  for (const unsigned int j : GeometryInfo<dim>::vertex_indices())
3077  out << cell->vertex(j)[i] << " ";
3078  }
3079  out << '\n';
3080  }
3081 
3082  // (v) write out cells to plot as quadplot objects
3083  out << "\n#"
3084  << "\n# List of cells to quadplot."
3085  << "\n# quadplot <vertex order> <id> <style>"
3086  << "\n#" << '\n';
3087  for (unsigned int i = 0; i < tria.n_active_cells(); ++i)
3088  {
3089  out << "\nquadplot f ";
3090  for (unsigned int j = 0; j < dim; ++j)
3091  out << axes[j] << i << " ";
3092  out << "\'k#\'";
3093  }
3094  out << '\n';
3095 
3096  // (vi) write footer
3097  out << "\n#"
3098  << "\n#"
3099  << "\n#" << '\n';
3100 
3101  // make sure everything now gets to the output stream
3102  out.flush();
3103  AssertThrow(out.fail() == false, ExcIO());
3104 }
3105 
3106 
3107 
3108 namespace
3109 {
3116  template <int dim, int spacedim, typename ITERATOR, typename END>
3117  void
3118  generate_triangulation_patches(
3119  std::vector<DataOutBase::Patch<dim, spacedim>> &patches,
3120  ITERATOR cell,
3121  END end)
3122  {
3123  // convert each of the active cells into a patch
3124  for (; cell != end; ++cell)
3125  {
3127  patch.reference_cell = cell->reference_cell();
3128  patch.n_subdivisions = 1;
3129  patch.data.reinit(5, cell->n_vertices());
3130 
3131  for (const unsigned int v : cell->vertex_indices())
3132  {
3133  patch.vertices[v] = cell->vertex(v);
3134  patch.data(0, v) = cell->level();
3135  patch.data(1, v) =
3136  static_cast<std::make_signed<types::manifold_id>::type>(
3137  cell->manifold_id());
3138  patch.data(2, v) =
3139  static_cast<std::make_signed<types::material_id>::type>(
3140  cell->material_id());
3141  if (cell->is_active())
3142  patch.data(3, v) =
3143  static_cast<std::make_signed<types::subdomain_id>::type>(
3144  cell->subdomain_id());
3145  else
3146  patch.data(3, v) = -1;
3147  patch.data(4, v) =
3148  static_cast<std::make_signed<types::subdomain_id>::type>(
3149  cell->level_subdomain_id());
3150  }
3151  patches.push_back(patch);
3152  }
3153  }
3154 
3155 
3156 
3157  std::vector<std::string>
3158  triangulation_patch_data_names()
3159  {
3160  std::vector<std::string> v(5);
3161  v[0] = "level";
3162  v[1] = "manifold";
3163  v[2] = "material";
3164  v[3] = "subdomain";
3165  v[4] = "level_subdomain";
3166  return v;
3167  }
3168 
3172  std::vector<typename Triangulation<3, 3>::active_line_iterator>
3173  get_boundary_edge_iterators(const Triangulation<3, 3> &tria)
3174  {
3175  std::vector<typename Triangulation<3, 3>::active_line_iterator> res;
3176 
3177  std::vector<bool> flags;
3178  tria.save_user_flags_line(flags);
3179  const_cast<Triangulation<3, 3> &>(tria).clear_user_flags_line();
3180 
3181  for (auto face : tria.active_face_iterators())
3182  for (const auto l : face->line_indices())
3183  {
3184  const auto line = face->line(l);
3185  if (line->user_flag_set() || line->has_children())
3186  continue;
3187  else
3188  line->set_user_flag();
3189  if (line->at_boundary())
3190  res.emplace_back(line);
3191  }
3192  const_cast<Triangulation<3, 3> &>(tria).load_user_flags_line(flags);
3193  return res;
3194  }
3195 
3196 
3197 
3201  template <int dim, int spacedim>
3202  std::vector<typename Triangulation<dim, spacedim>::active_line_iterator>
3203  get_boundary_edge_iterators(const Triangulation<dim, spacedim> &)
3204  {
3205  return {};
3206  }
3207 
3208 
3209 
3214  std::vector<typename Triangulation<3, 3>::active_line_iterator>
3215  get_relevant_edge_iterators(const Triangulation<3, 3> &tria)
3216  {
3217  std::vector<typename Triangulation<3, 3>::active_line_iterator> res;
3218 
3219  std::vector<bool> flags;
3220  tria.save_user_flags_line(flags);
3221  const_cast<Triangulation<3, 3> &>(tria).clear_user_flags_line();
3222 
3223  for (auto face : tria.active_face_iterators())
3224  for (const auto l : face->line_indices())
3225  {
3226  const auto line = face->line(l);
3227  if (line->user_flag_set() || line->has_children())
3228  continue;
3229  else
3230  line->set_user_flag();
3231  if (line->manifold_id() != numbers::flat_manifold_id ||
3232  (line->boundary_id() != 0 &&
3233  line->boundary_id() != numbers::invalid_boundary_id))
3234  res.emplace_back(line);
3235  }
3236  const_cast<Triangulation<3, 3> &>(tria).load_user_flags_line(flags);
3237  return res;
3238  }
3239 
3240 
3244  template <int dim, int spacedim>
3245  std::vector<typename Triangulation<dim, spacedim>::active_line_iterator>
3246  get_relevant_edge_iterators(const Triangulation<dim, spacedim> &)
3247  {
3248  return {};
3249  }
3250 
3251 
3252 
3256  template <int dim, int spacedim>
3257  std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3258  get_boundary_face_iterators(const Triangulation<dim, spacedim> &tria)
3259  {
3260  std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3261  res;
3262  if (dim == 1)
3263  return res;
3264  for (auto face : tria.active_face_iterators())
3265  {
3266  if (face->boundary_id() != numbers::invalid_boundary_id)
3267  res.push_back(face);
3268  }
3269  return res;
3270  }
3271 
3272 
3273 
3278  template <int dim, int spacedim>
3279  std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3280  get_relevant_face_iterators(const Triangulation<dim, spacedim> &tria)
3281  {
3282  std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3283  res;
3284  if (dim == 1)
3285  return res;
3286  for (auto face : tria.active_face_iterators())
3287  {
3288  if (face->manifold_id() != numbers::flat_manifold_id ||
3289  (face->boundary_id() != 0 &&
3290  face->boundary_id() != numbers::invalid_boundary_id))
3291  res.push_back(face);
3292  }
3293  return res;
3294  }
3295 } // namespace
3296 
3297 
3298 
3299 template <int dim, int spacedim>
3300 void
3302  std::ostream & out) const
3303 {
3304  AssertThrow(out.fail() == false, ExcIO());
3305 
3306  // get the positions of the vertices
3307  const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
3308 
3309  const auto n_vertices = vertices.size();
3310 
3311  out << "# vtk DataFile Version 3.0\n"
3312  << "Triangulation generated with deal.II\n"
3313  << "ASCII\n"
3314  << "DATASET UNSTRUCTURED_GRID\n"
3315  << "POINTS " << n_vertices << " double\n";
3316 
3317  // actually write the vertices.
3318  for (const auto &v : vertices)
3319  {
3320  out << v;
3321  for (unsigned int d = spacedim + 1; d <= 3; ++d)
3322  out << " 0"; // fill with zeroes
3323  out << '\n';
3324  }
3325 
3326  const auto faces = vtk_flags.output_only_relevant ?
3327  get_relevant_face_iterators(tria) :
3328  get_boundary_face_iterators(tria);
3329  const auto edges = vtk_flags.output_only_relevant ?
3330  get_relevant_edge_iterators(tria) :
3331  get_boundary_edge_iterators(tria);
3332 
3333  AssertThrow(
3334  vtk_flags.output_cells || (dim >= 2 && vtk_flags.output_faces) ||
3335  (dim >= 3 && vtk_flags.output_edges),
3336  ExcMessage(
3337  "At least one of the flags (output_cells, output_faces, output_edges) has to be enabled!"));
3338 
3339  // Write cells preamble
3340  const int n_cells = (vtk_flags.output_cells ? tria.n_active_cells() : 0) +
3341  (vtk_flags.output_faces ? faces.size() : 0) +
3342  (vtk_flags.output_edges ? edges.size() : 0);
3343 
3344  // VTK now expects a number telling the total storage requirement to read all
3345  // cell connectivity information. The connectivity information is read cell by
3346  // cell, first specifying how many vertices are required to describe the cell,
3347  // and then specifying the index of every vertex. This means that for every
3348  // deal.II object type, we always need n_vertices + 1 integer per cell.
3349  // Compute the total number here.
3350  int cells_size = 0;
3351 
3352  if (vtk_flags.output_cells)
3353  for (const auto &cell : tria.active_cell_iterators())
3354  cells_size += cell->n_vertices() + 1;
3355 
3356  if (vtk_flags.output_faces)
3357  for (const auto &face : faces)
3358  cells_size += face->n_vertices() + 1;
3359 
3360  if (vtk_flags.output_edges)
3361  for (const auto &edge : edges)
3362  cells_size += edge->n_vertices() + 1;
3363 
3364  AssertThrow(cells_size > 0, ExcMessage("No cells given to be output!"));
3365 
3366  out << "\nCELLS " << n_cells << ' ' << cells_size << '\n';
3367  /*
3368  * VTK cells:
3369  *
3370  * 1 VTK_VERTEX
3371  * 3 VTK_LINE
3372  * 5 VTK_TRIANGLE
3373  * 9 VTK_QUAD
3374  * 10 VTK_TETRA
3375  * 14 VTK_PYRAMID
3376  * 13 VTK_WEDGE
3377  * 12 VTK_HEXAHEDRON
3378  *
3379  * see also: https://vtk.org/wp-content/uploads/2015/04/file-formats.pdf
3380  */
3381  static const std::array<int, 8> deal_to_vtk_cell_type = {
3382  {1, 3, 5, 9, 10, 14, 13, 12}};
3383  static const std::array<unsigned int, 8> vtk_to_deal_hypercube = {
3384  {0, 1, 3, 2, 4, 5, 7, 6}};
3385 
3386  // write cells.
3387  if (vtk_flags.output_cells)
3388  for (const auto &cell : tria.active_cell_iterators())
3389  {
3390  out << cell->n_vertices();
3391  for (const unsigned int i : cell->vertex_indices())
3392  {
3393  out << ' ';
3394  const auto reference_cell = cell->reference_cell();
3395 
3400  out << cell->vertex_index(vtk_to_deal_hypercube[i]);
3401  else if ((reference_cell == ReferenceCells::Triangle) ||
3404  out << cell->vertex_index(i);
3406  {
3407  static const std::array<unsigned int, 5> permutation_table{
3408  {0, 1, 3, 2, 4}};
3409  out << cell->vertex_index(permutation_table[i]);
3410  }
3411  else
3412  Assert(false, ExcNotImplemented());
3413  }
3414  out << '\n';
3415  }
3416  if (vtk_flags.output_faces)
3417  for (const auto &face : faces)
3418  {
3419  out << face->n_vertices();
3420  for (const unsigned int i : face->vertex_indices())
3421  {
3422  out << ' '
3423  << face->vertex_index(GeometryInfo<dim>::vertices_per_face ==
3424  face->n_vertices() ?
3425  vtk_to_deal_hypercube[i] :
3426  i);
3427  }
3428  out << '\n';
3429  }
3430  if (vtk_flags.output_edges)
3431  for (const auto &edge : edges)
3432  {
3433  out << 2;
3434  for (const unsigned int i : edge->vertex_indices())
3435  out << ' ' << edge->vertex_index(i);
3436  out << '\n';
3437  }
3438 
3439  // write cell types
3440  out << "\nCELL_TYPES " << n_cells << '\n';
3441  if (vtk_flags.output_cells)
3442  {
3443  for (const auto &cell : tria.active_cell_iterators())
3444  out << deal_to_vtk_cell_type[static_cast<int>(cell->reference_cell())]
3445  << ' ';
3446  out << '\n';
3447  }
3448  if (vtk_flags.output_faces)
3449  {
3450  for (const auto &face : faces)
3451  out << deal_to_vtk_cell_type[static_cast<int>(face->reference_cell())]
3452  << ' ';
3453  out << '\n';
3454  }
3455  if (vtk_flags.output_edges)
3456  {
3457  for (const auto &edge : edges)
3458  out << deal_to_vtk_cell_type[static_cast<int>(edge->reference_cell())]
3459  << ' ';
3460  }
3461  out << "\n\nCELL_DATA " << n_cells << '\n'
3462  << "SCALARS MaterialID int 1\n"
3463  << "LOOKUP_TABLE default\n";
3464 
3465  // Now material id and boundary id
3466  if (vtk_flags.output_cells)
3467  {
3468  for (const auto &cell : tria.active_cell_iterators())
3469  {
3470  out << static_cast<std::make_signed<types::material_id>::type>(
3471  cell->material_id())
3472  << ' ';
3473  }
3474  out << '\n';
3475  }
3476  if (vtk_flags.output_faces)
3477  {
3478  for (const auto &face : faces)
3479  {
3480  out << static_cast<std::make_signed<types::boundary_id>::type>(
3481  face->boundary_id())
3482  << ' ';
3483  }
3484  out << '\n';
3485  }
3486  if (vtk_flags.output_edges)
3487  {
3488  for (const auto &edge : edges)
3489  {
3490  out << static_cast<std::make_signed<types::boundary_id>::type>(
3491  edge->boundary_id())
3492  << ' ';
3493  }
3494  }
3495 
3496  out << "\n\nSCALARS ManifoldID int 1\n"
3497  << "LOOKUP_TABLE default\n";
3498 
3499  // Now manifold id
3500  if (vtk_flags.output_cells)
3501  {
3502  for (const auto &cell : tria.active_cell_iterators())
3503  {
3504  out << static_cast<std::make_signed<types::manifold_id>::type>(
3505  cell->manifold_id())
3506  << ' ';
3507  }
3508  out << '\n';
3509  }
3510  if (vtk_flags.output_faces)
3511  {
3512  for (const auto &face : faces)
3513  {
3514  out << static_cast<std::make_signed<types::manifold_id>::type>(
3515  face->manifold_id())
3516  << ' ';
3517  }
3518  out << '\n';
3519  }
3520  if (vtk_flags.output_edges)
3521  {
3522  for (const auto &edge : edges)
3523  {
3524  out << static_cast<std::make_signed<types::manifold_id>::type>(
3525  edge->manifold_id())
3526  << ' ';
3527  }
3528  out << '\n';
3529  }
3530 
3531  out.flush();
3532 
3533  AssertThrow(out.fail() == false, ExcIO());
3534 }
3535 
3536 
3537 
3538 template <int dim, int spacedim>
3539 void
3541  std::ostream & out) const
3542 {
3543  AssertThrow(out.fail() == false, ExcIO());
3544 
3545  // convert the cells of the triangulation into a set of patches
3546  // and then have them output. since there is no data attached to
3547  // the geometry, we also do not have to provide any names, identifying
3548  // information, etc.
3549  std::vector<DataOutBase::Patch<dim, spacedim>> patches;
3550  patches.reserve(tria.n_active_cells());
3551  generate_triangulation_patches(patches, tria.begin_active(), tria.end());
3552 
3555  patches,
3556  triangulation_patch_data_names(),
3557  std::vector<
3558  std::tuple<unsigned int,
3559  unsigned int,
3560  std::string,
3562  vtu_flags,
3563  out);
3565  {
3566  out << " </UnstructuredGrid>\n";
3567  out << "<dealiiData encoding=\"base64\">";
3568  std::stringstream outstring;
3569  boost::archive::binary_oarchive ia(outstring);
3570  tria.save(ia, 0);
3571  const auto compressed = Utilities::compress(outstring.str());
3572  out << Utilities::encode_base64({compressed.begin(), compressed.end()});
3573  out << "\n</dealiiData>\n";
3574  out << "</VTKFile>\n";
3575  }
3576  else
3578 
3579  out << std::flush;
3580  AssertThrow(out.fail() == false, ExcIO());
3581 }
3582 
3583 
3584 
3585 template <int dim, int spacedim>
3586 void
3589  const std::string & filename_without_extension,
3590  const bool view_levels,
3591  const bool include_artificial) const
3592 {
3593  std::vector<DataOutBase::Patch<dim, spacedim>> patches;
3594  const unsigned int n_datasets = 4;
3595  std::vector<std::string> data_names;
3596  data_names.emplace_back("level");
3597  data_names.emplace_back("subdomain");
3598  data_names.emplace_back("level_subdomain");
3599  data_names.emplace_back("proc_writing");
3600 
3601  const auto reference_cells = tria.get_reference_cells();
3602 
3603  AssertDimension(reference_cells.size(), 1);
3604 
3605  const auto &reference_cell = reference_cells[0];
3606 
3607  const unsigned int n_q_points = reference_cell.n_vertices();
3608 
3609  for (const auto &cell : tria.cell_iterators())
3610  {
3611  if (!view_levels)
3612  {
3613  if (cell->has_children())
3614  continue;
3615  if (!include_artificial &&
3616  cell->subdomain_id() == numbers::artificial_subdomain_id)
3617  continue;
3618  }
3619  else if (!include_artificial)
3620  {
3621  if (cell->has_children() &&
3622  cell->level_subdomain_id() == numbers::artificial_subdomain_id)
3623  continue;
3624  else if (cell->is_active() &&
3625  cell->level_subdomain_id() ==
3627  cell->subdomain_id() == numbers::artificial_subdomain_id)
3628  continue;
3629  }
3630 
3632  patch.data.reinit(n_datasets, n_q_points);
3633  patch.points_are_available = false;
3635 
3636  for (unsigned int vertex = 0; vertex < n_q_points; ++vertex)
3637  {
3638  patch.vertices[vertex] = cell->vertex(vertex);
3639  patch.data(0, vertex) = cell->level();
3640  if (cell->is_active())
3641  patch.data(1, vertex) = static_cast<double>(
3642  static_cast<std::make_signed<types::subdomain_id>::type>(
3643  cell->subdomain_id()));
3644  else
3645  patch.data(1, vertex) = -1.0;
3646  patch.data(2, vertex) = static_cast<double>(
3647  static_cast<std::make_signed<types::subdomain_id>::type>(
3648  cell->level_subdomain_id()));
3649  patch.data(3, vertex) = tria.locally_owned_subdomain();
3650  }
3651 
3652  for (auto f : reference_cell.face_indices())
3654  patches.push_back(patch);
3655  }
3656 
3657  // only create .pvtu file if running in parallel
3658  // if not, just create a .vtu file with no reference
3659  // to the processor number
3660  std::string new_file = filename_without_extension + ".vtu";
3662  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(&tria))
3663  {
3664  new_file = filename_without_extension + ".proc" +
3665  Utilities::int_to_string(tr->locally_owned_subdomain(), 4) +
3666  ".vtu";
3667 
3668  // create .pvtu record
3669  if (tr->locally_owned_subdomain() == 0)
3670  {
3671  std::vector<std::string> filenames;
3672 
3673  // .pvtu needs to reference the files without a relative path because
3674  // it will be written in the same directory. For this, remove any
3675  // paths from filename.
3676  std::size_t pos = filename_without_extension.find_last_of('/');
3677  if (pos == std::string::npos)
3678  pos = 0;
3679  else
3680  pos += 1;
3681  const unsigned int n_procs =
3682  Utilities::MPI::n_mpi_processes(tr->get_communicator());
3683  for (unsigned int i = 0; i < n_procs; ++i)
3684  filenames.push_back(filename_without_extension.substr(pos) +
3685  ".proc" + Utilities::int_to_string(i, 4) +
3686  ".vtu");
3687 
3688  const std::string pvtu_filename =
3689  (filename_without_extension + ".pvtu");
3690  std::ofstream pvtu_output(pvtu_filename.c_str());
3691 
3692  DataOut<dim, spacedim> data_out;
3693  data_out.attach_triangulation(*tr);
3694 
3695  // We need a dummy vector with the names of the data values in the
3696  // .vtu files in order that the .pvtu contains reference these values
3697  Vector<float> dummy_vector(tr->n_active_cells());
3698  data_out.add_data_vector(dummy_vector, "level");
3699  data_out.add_data_vector(dummy_vector, "subdomain");
3700  data_out.add_data_vector(dummy_vector, "level_subdomain");
3701  data_out.add_data_vector(dummy_vector, "proc_writing");
3702 
3703  data_out.build_patches();
3704 
3705  data_out.write_pvtu_record(pvtu_output, filenames);
3706  }
3707  }
3708 
3709  std::ofstream out(new_file.c_str());
3710  std::vector<
3711  std::tuple<unsigned int,
3712  unsigned int,
3713  std::string,
3715  vector_data_ranges;
3716  DataOutBase::VtkFlags flags;
3717  DataOutBase::write_vtu(patches, data_names, vector_data_ranges, flags, out);
3718 }
3719 
3720 
3721 
3722 unsigned int
3724 {
3725  return 0;
3726 }
3727 
3728 unsigned int
3730 {
3731  return 0;
3732 }
3733 
3734 
3735 unsigned int
3737 {
3738  return 0;
3739 }
3740 
3741 unsigned int
3743 {
3744  return 0;
3745 }
3746 
3747 unsigned int
3749 {
3750  return 0;
3751 }
3752 
3753 unsigned int
3755 {
3756  return 0;
3757 }
3758 
3759 unsigned int
3761 {
3762  return 0;
3763 }
3764 
3765 unsigned int
3767 {
3768  return 0;
3769 }
3770 
3771 
3772 
3773 template <int dim, int spacedim>
3774 unsigned int
3776 {
3778  unsigned int n_faces = 0;
3779 
3780  for (const auto &face : tria.active_face_iterators())
3781  if ((face->at_boundary()) && (face->boundary_id() != 0))
3782  n_faces++;
3783 
3784  return n_faces;
3785 }
3786 
3787 
3788 
3789 template <int dim, int spacedim>
3790 unsigned int
3792 {
3793  // save the user flags for lines so
3794  // we can use these flags to track
3795  // which ones we've already counted
3796  std::vector<bool> line_flags;
3797  const_cast<::Triangulation<dim, spacedim> &>(tria).save_user_flags_line(
3798  line_flags);
3799  const_cast<::Triangulation<dim, spacedim> &>(tria)
3800  .clear_user_flags_line();
3801 
3802  unsigned int n_lines = 0;
3803 
3804  for (const auto &cell : tria.active_cell_iterators())
3805  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
3806  if (cell->line(l)->at_boundary() && (cell->line(l)->boundary_id() != 0) &&
3807  (cell->line(l)->user_flag_set() == false))
3808  {
3809  ++n_lines;
3810  cell->line(l)->set_user_flag();
3811  }
3812 
3813  // at the end, restore the user
3814  // flags for the lines
3815  const_cast<::Triangulation<dim, spacedim> &>(tria).load_user_flags_line(
3816  line_flags);
3817 
3818  return n_lines;
3819 }
3820 
3821 
3822 
3823 unsigned int
3825  const unsigned int next_element_index,
3826  std::ostream &) const
3827 {
3828  return next_element_index;
3829 }
3830 
3831 
3832 unsigned int
3834  const unsigned int next_element_index,
3835  std::ostream &) const
3836 {
3837  return next_element_index;
3838 }
3839 
3840 unsigned int
3842  const unsigned int next_element_index,
3843  std::ostream &) const
3844 {
3845  return next_element_index;
3846 }
3847 
3848 
3849 unsigned int
3851  const unsigned int next_element_index,
3852  std::ostream &) const
3853 {
3854  return next_element_index;
3855 }
3856 
3857 unsigned int
3859  const unsigned int next_element_index,
3860  std::ostream &) const
3861 {
3862  return next_element_index;
3863 }
3864 
3865 
3866 unsigned int
3868  const unsigned int next_element_index,
3869  std::ostream &) const
3870 {
3871  return next_element_index;
3872 }
3873 
3874 
3875 unsigned int
3877  const unsigned int next_element_index,
3878  std::ostream &) const
3879 {
3880  return next_element_index;
3881 }
3882 
3883 unsigned int
3885  const unsigned int next_element_index,
3886  std::ostream &) const
3887 {
3888  return next_element_index;
3889 }
3890 
3891 
3892 
3893 template <int dim, int spacedim>
3894 unsigned int
3896  const unsigned int next_element_index,
3897  std::ostream & out) const
3898 {
3899  unsigned int current_element_index = next_element_index;
3900 
3901  for (const auto &face : tria.active_face_iterators())
3902  if (face->at_boundary() && (face->boundary_id() != 0))
3903  {
3904  out << current_element_index << ' '
3905  << face->reference_cell().gmsh_element_type() << ' ';
3906  out << static_cast<unsigned int>(face->boundary_id()) << ' '
3907  << static_cast<unsigned int>(face->boundary_id()) << ' '
3908  << face->n_vertices();
3909  // note: vertex numbers are 1-base
3910  for (unsigned int vertex : face->vertex_indices())
3911  {
3912  if (face->reference_cell() == ReferenceCells::Quadrilateral)
3913  out << ' '
3914  << face->vertex_index(
3916  1;
3917  else if ((face->reference_cell() == ReferenceCells::Triangle) ||
3918  (face->reference_cell() == ReferenceCells::Line))
3919  out << ' ' << face->vertex_index(vertex) + 1;
3920  else
3921  Assert(false, ExcInternalError());
3922  }
3923  out << '\n';
3924 
3925  ++current_element_index;
3926  }
3927  return current_element_index;
3928 }
3929 
3930 
3931 
3932 template <int dim, int spacedim>
3933 unsigned int
3935  const unsigned int next_element_index,
3936  std::ostream & out) const
3937 {
3938  unsigned int current_element_index = next_element_index;
3939  // save the user flags for lines so
3940  // we can use these flags to track
3941  // which ones we've already taken
3942  // care of
3943  std::vector<bool> line_flags;
3944  const_cast<::Triangulation<dim, spacedim> &>(tria).save_user_flags_line(
3945  line_flags);
3946  const_cast<::Triangulation<dim, spacedim> &>(tria)
3947  .clear_user_flags_line();
3948 
3949  for (const auto &cell : tria.active_cell_iterators())
3950  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
3951  if (cell->line(l)->at_boundary() && (cell->line(l)->boundary_id() != 0) &&
3952  (cell->line(l)->user_flag_set() == false))
3953  {
3954  out << next_element_index << ' '
3956  out << static_cast<unsigned int>(cell->line(l)->boundary_id()) << ' '
3957  << static_cast<unsigned int>(cell->line(l)->boundary_id())
3958  << " 2 "; // two vertex indices to follow
3959  // note: vertex numbers are 1-base
3960  for (unsigned int vertex = 0; vertex < 2; ++vertex)
3961  out << ' '
3962  << cell->line(l)->vertex_index(
3964  1;
3965  out << '\n';
3966 
3967  // move on to the next line
3968  // but mark the current one
3969  // as taken care of
3970  ++current_element_index;
3971  cell->line(l)->set_user_flag();
3972  }
3973 
3974  // at the end, restore the user
3975  // flags for the lines
3976  const_cast<::Triangulation<dim, spacedim> &>(tria).load_user_flags_line(
3977  line_flags);
3978 
3979  return current_element_index;
3980 }
3981 
3982 
3983 
3984 unsigned int
3986  const unsigned int next_element_index,
3987  std::ostream &) const
3988 {
3989  return next_element_index;
3990 }
3991 
3992 unsigned int
3994  const unsigned int next_element_index,
3995  std::ostream &) const
3996 {
3997  return next_element_index;
3998 }
3999 
4000 unsigned int
4002  const unsigned int next_element_index,
4003  std::ostream &) const
4004 {
4005  return next_element_index;
4006 }
4007 
4008 unsigned int
4010  const unsigned int next_element_index,
4011  std::ostream &) const
4012 {
4013  return next_element_index;
4014 }
4015 
4016 unsigned int
4018  const unsigned int next_element_index,
4019  std::ostream &) const
4020 {
4021  return next_element_index;
4022 }
4023 
4024 
4025 unsigned int
4027  const unsigned int next_element_index,
4028  std::ostream &) const
4029 {
4030  return next_element_index;
4031 }
4032 
4033 
4034 unsigned int
4036  const unsigned int next_element_index,
4037  std::ostream &) const
4038 {
4039  return next_element_index;
4040 }
4041 
4042 unsigned int
4044  const unsigned int next_element_index,
4045  std::ostream &) const
4046 {
4047  return next_element_index;
4048 }
4049 
4050 
4051 
4052 template <int dim, int spacedim>
4053 unsigned int
4055  const unsigned int next_element_index,
4056  std::ostream & out) const
4057 {
4058  unsigned int current_element_index = next_element_index;
4060 
4061  for (const auto &face : tria.active_face_iterators())
4062  if (face->at_boundary() && (face->boundary_id() != 0))
4063  {
4064  out << current_element_index << " "
4065  << static_cast<unsigned int>(face->boundary_id()) << " ";
4066  switch (dim)
4067  {
4068  case 2:
4069  out << "line ";
4070  break;
4071  case 3:
4072  out << "quad ";
4073  break;
4074  default:
4075  Assert(false, ExcNotImplemented());
4076  }
4077  // note: vertex numbers are 1-base
4078  for (unsigned int vertex = 0;
4079  vertex < GeometryInfo<dim>::vertices_per_face;
4080  ++vertex)
4081  out << face->vertex_index(
4083  1
4084  << ' ';
4085  out << '\n';
4086 
4087  ++current_element_index;
4088  }
4089  return current_element_index;
4090 }
4091 
4092 
4093 
4094 template <int dim, int spacedim>
4095 unsigned int
4097  const unsigned int next_element_index,
4098  std::ostream & out) const
4099 {
4100  unsigned int current_element_index = next_element_index;
4101  // save the user flags for lines so
4102  // we can use these flags to track
4103  // which ones we've already taken
4104  // care of
4105  std::vector<bool> line_flags;
4106  const_cast<::Triangulation<dim, spacedim> &>(tria).save_user_flags_line(
4107  line_flags);
4108  const_cast<::Triangulation<dim, spacedim> &>(tria)
4109  .clear_user_flags_line();
4110 
4111  for (const auto &cell : tria.active_cell_iterators())
4112  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
4113  if (cell->line(l)->at_boundary() && (cell->line(l)->boundary_id() != 0) &&
4114  (cell->line(l)->user_flag_set() == false))
4115  {
4116  out << current_element_index << " "
4117  << static_cast<unsigned int>(cell->line(l)->boundary_id())
4118  << " line ";
4119  // note: vertex numbers in ucd format are 1-base
4120  for (unsigned int vertex = 0; vertex < 2; ++vertex)
4121  out << cell->line(l)->vertex_index(
4123  1
4124  << ' ';
4125  out << '\n';
4126 
4127  // move on to the next line
4128  // but mark the current one
4129  // as taken care of
4130  ++current_element_index;
4131  cell->line(l)->set_user_flag();
4132  }
4133 
4134  // at the end, restore the user
4135  // flags for the lines
4136  const_cast<::Triangulation<dim, spacedim> &>(tria).load_user_flags_line(
4137  line_flags);
4138  return current_element_index;
4139 }
4140 
4141 
4142 
4143 namespace internal
4144 {
4145  namespace
4146  {
4155  template <int spacedim>
4156  void
4157  remove_colinear_points(std::vector<Point<spacedim>> &points)
4158  {
4159  while (points.size() > 2)
4160  {
4161  Tensor<1, spacedim> first_difference = points[1] - points[0];
4162  first_difference /= first_difference.norm();
4163  Tensor<1, spacedim> second_difference = points[2] - points[1];
4164  second_difference /= second_difference.norm();
4165  // If the three points are colinear then remove the middle one.
4166  if ((first_difference - second_difference).norm() < 1e-10)
4167  points.erase(points.begin() + 1);
4168  else
4169  break;
4170  }
4171  }
4172 
4173 
4174 
4175  template <int spacedim>
4176  void
4177  write_gnuplot(const ::Triangulation<1, spacedim> &tria,
4178  std::ostream & out,
4179  const Mapping<1, spacedim> *,
4180  const GridOutFlags::Gnuplot &gnuplot_flags)
4181  {
4182  AssertThrow(out.fail() == false, ExcIO());
4183 
4184  for (const auto &cell : tria.active_cell_iterators())
4185  {
4186  if (gnuplot_flags.write_cell_numbers)
4187  out << "# cell " << cell << '\n';
4188 
4189  out << cell->vertex(0) << ' ' << cell->level() << ' '
4190  << cell->material_id() << '\n'
4191  << cell->vertex(1) << ' ' << cell->level() << ' '
4192  << cell->material_id() << '\n'
4193  << "\n\n";
4194  }
4195 
4196  // make sure everything now gets to
4197  // disk
4198  out.flush();
4199 
4200  AssertThrow(out.fail() == false, ExcIO());
4201  }
4202 
4203 
4204 
4205  template <int spacedim>
4206  void
4207  write_gnuplot(const ::Triangulation<2, spacedim> &tria,
4208  std::ostream & out,
4209  const Mapping<2, spacedim> * mapping,
4210  const GridOutFlags::Gnuplot & gnuplot_flags)
4211  {
4212  AssertThrow(out.fail() == false, ExcIO());
4213 
4214  const int dim = 2;
4215 
4216  const unsigned int n_additional_points =
4217  gnuplot_flags.n_extra_curved_line_points;
4218  const unsigned int n_points = 2 + n_additional_points;
4219 
4220  // If we need to plot curved lines then generate a quadrature formula to
4221  // place points via the mapping
4222  Quadrature<dim> q_projector;
4223  std::vector<Point<dim - 1>> boundary_points;
4224  if (mapping != nullptr)
4225  {
4226  boundary_points.resize(n_points);
4227  boundary_points[0][0] = 0;
4228  boundary_points[n_points - 1][0] = 1;
4229  for (unsigned int i = 1; i < n_points - 1; ++i)
4230  boundary_points[i](0) = 1. * i / (n_points - 1);
4231 
4232  std::vector<double> dummy_weights(n_points, 1. / n_points);
4233  Quadrature<dim - 1> quadrature(boundary_points, dummy_weights);
4234 
4236  ::ReferenceCells::Quadrilateral, quadrature);
4237  }
4238 
4239  for (const auto &cell : tria.active_cell_iterators())
4240  {
4241  if (gnuplot_flags.write_cell_numbers)
4242  out << "# cell " << cell << '\n';
4243 
4244  if (mapping == nullptr ||
4245  (dim == spacedim ?
4246  (!cell->at_boundary() && !gnuplot_flags.curved_inner_cells) :
4247  // ignore checking for boundary or interior cells in the codim
4248  // 1 case: 'or false' is a no-op
4249  false))
4250  {
4251  // write out the four sides of this cell by putting the four
4252  // points (+ the initial point again) in a row and lifting the
4253  // drawing pencil at the end
4254  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
4255  out << cell->vertex(GeometryInfo<dim>::ucd_to_deal[i]) << ' '
4256  << cell->level() << ' ' << cell->material_id() << '\n';
4257  out << cell->vertex(0) << ' ' << cell->level() << ' '
4258  << cell->material_id() << '\n'
4259  << '\n' // double new line for gnuplot 3d plots
4260  << '\n';
4261  }
4262  else
4263  // cell is at boundary and we are to treat curved boundaries. so
4264  // loop over all faces and draw them as small pieces of lines
4265  {
4266  for (const unsigned int face_no :
4268  {
4269  const typename ::Triangulation<dim,
4270  spacedim>::face_iterator
4271  face = cell->face(face_no);
4272  if (dim != spacedim || face->at_boundary() ||
4273  gnuplot_flags.curved_inner_cells)
4274  {
4275  // Save the points on each face to a vector and then try
4276  // to remove colinear points that won't show up in the
4277  // generated plot.
4278  std::vector<Point<spacedim>> line_points;
4279  // compute offset of quadrature points within set of
4280  // projected points
4281  const unsigned int offset = face_no * n_points;
4282  for (unsigned int i = 0; i < n_points; ++i)
4283  line_points.push_back(
4284  mapping->transform_unit_to_real_cell(
4285  cell, q_projector.point(offset + i)));
4286  internal::remove_colinear_points(line_points);
4287 
4288  for (const Point<spacedim> &point : line_points)
4289  out << point << ' ' << cell->level() << ' '
4290  << cell->material_id() << '\n';
4291 
4292  out << '\n' << '\n';
4293  }
4294  else
4295  {
4296  // if, however, the face is not at the boundary and we
4297  // don't want to curve anything, then draw it as usual
4298  out << face->vertex(0) << ' ' << cell->level() << ' '
4299  << cell->material_id() << '\n'
4300  << face->vertex(1) << ' ' << cell->level() << ' '
4301  << cell->material_id() << '\n'
4302  << '\n'
4303  << '\n';
4304  }
4305  }
4306  }
4307  }
4308 
4309  // make sure everything now gets to disk
4310  out.flush();
4311 
4312  AssertThrow(out.fail() == false, ExcIO());
4313  }
4314 
4315 
4316 
4317  template <int spacedim>
4318  void
4319  write_gnuplot(const ::Triangulation<3, spacedim> &tria,
4320  std::ostream & out,
4321  const Mapping<3, spacedim> * mapping,
4322  const GridOutFlags::Gnuplot & gnuplot_flags)
4323  {
4324  AssertThrow(out.fail() == false, ExcIO());
4325 
4326  const int dim = 3;
4327 
4328  const unsigned int n_additional_points =
4329  gnuplot_flags.n_extra_curved_line_points;
4330  const unsigned int n_points = 2 + n_additional_points;
4331 
4332  // If we need to plot curved lines then generate a quadrature formula to
4333  // place points via the mapping
4334  std::unique_ptr<Quadrature<dim>> q_projector;
4335  std::vector<Point<1>> boundary_points;
4336  if (mapping != nullptr)
4337  {
4338  boundary_points.resize(n_points);
4339  boundary_points[0][0] = 0;
4340  boundary_points[n_points - 1][0] = 1;
4341  for (unsigned int i = 1; i < n_points - 1; ++i)
4342  boundary_points[i](0) = 1. * i / (n_points - 1);
4343 
4344  std::vector<double> dummy_weights(n_points, 1. / n_points);
4345  Quadrature<1> quadrature1d(boundary_points, dummy_weights);
4346 
4347  // tensor product of points, only one copy
4348  QIterated<dim - 1> quadrature(quadrature1d, 1);
4349  q_projector = std::make_unique<Quadrature<dim>>(
4351  }
4352 
4353  for (const auto &cell : tria.active_cell_iterators())
4354  {
4355  if (gnuplot_flags.write_cell_numbers)
4356  out << "# cell " << cell << '\n';
4357 
4358  if (mapping == nullptr || n_points == 2 ||
4359  (!cell->has_boundary_lines() &&
4360  !gnuplot_flags.curved_inner_cells))
4361  {
4362  if (cell->reference_cell() == ReferenceCells::Hexahedron)
4363  {
4364  // front face
4365  out << cell->vertex(0) << ' ' << cell->level() << ' '
4366  << cell->material_id() << '\n'
4367  << cell->vertex(1) << ' ' << cell->level() << ' '
4368  << cell->material_id() << '\n'
4369  << cell->vertex(5) << ' ' << cell->level() << ' '
4370  << cell->material_id() << '\n'
4371  << cell->vertex(4) << ' ' << cell->level() << ' '
4372  << cell->material_id() << '\n'
4373  << cell->vertex(0) << ' ' << cell->level() << ' '
4374  << cell->material_id() << '\n'
4375  << '\n';
4376  // back face
4377  out << cell->vertex(2) << ' ' << cell->level() << ' '
4378  << cell->material_id() << '\n'
4379  << cell->vertex(3) << ' ' << cell->level() << ' '
4380  << cell->material_id() << '\n'
4381  << cell->vertex(7) << ' ' << cell->level() << ' '
4382  << cell->material_id() << '\n'
4383  << cell->vertex(6) << ' ' << cell->level() << ' '
4384  << cell->material_id() << '\n'
4385  << cell->vertex(2) << ' ' << cell->level() << ' '
4386  << cell->material_id() << '\n'
4387  << '\n';
4388 
4389  // now for the four connecting lines
4390  out << cell->vertex(0) << ' ' << cell->level() << ' '
4391  << cell->material_id() << '\n'
4392  << cell->vertex(2) << ' ' << cell->level() << ' '
4393  << cell->material_id() << '\n'
4394  << '\n';
4395  out << cell->vertex(1) << ' ' << cell->level() << ' '
4396  << cell->material_id() << '\n'
4397  << cell->vertex(3) << ' ' << cell->level() << ' '
4398  << cell->material_id() << '\n'
4399  << '\n';
4400  out << cell->vertex(5) << ' ' << cell->level() << ' '
4401  << cell->material_id() << '\n'
4402  << cell->vertex(7) << ' ' << cell->level() << ' '
4403  << cell->material_id() << '\n'
4404  << '\n';
4405  out << cell->vertex(4) << ' ' << cell->level() << ' '
4406  << cell->material_id() << '\n'
4407  << cell->vertex(6) << ' ' << cell->level() << ' '
4408  << cell->material_id() << '\n'
4409  << '\n';
4410  }
4411  else if (cell->reference_cell() == ReferenceCells::Tetrahedron)
4412  {
4413  // Draw the tetrahedron as a two collections of lines.
4414  for (const unsigned int v : {0, 1, 2, 0, 3, 2})
4415  {
4416  out << cell->vertex(v) << ' ' << cell->level() << ' '
4417  << cell->material_id() << '\n';
4418  }
4419  out << '\n'; // end of first line
4420 
4421  for (const unsigned int v : {3, 1})
4422  {
4423  out << cell->vertex(v) << ' ' << cell->level() << ' '
4424  << cell->material_id() << '\n';
4425  }
4426  out << '\n'; // end of second line
4427  }
4428  else if (cell->reference_cell() == ReferenceCells::Wedge)
4429  {
4430  // Draw the wedge as a collection of three
4431  // lines. The first one wraps around the base,
4432  // goes up to the top, and wraps around that. The
4433  // second and third are just individual lines
4434  // going from base to top.
4435  for (const unsigned int v : {0, 1, 2, 0, 3, 4, 5, 3})
4436  {
4437  out << cell->vertex(v) << ' ' << cell->level() << ' '
4438  << cell->material_id() << '\n';
4439  }
4440  out << '\n'; // end of first line
4441 
4442  for (const unsigned int v : {1, 4})
4443  {
4444  out << cell->vertex(v) << ' ' << cell->level() << ' '
4445  << cell->material_id() << '\n';
4446  }
4447  out << '\n'; // end of second line
4448 
4449  for (const unsigned int v : {2, 5})
4450  {
4451  out << cell->vertex(v) << ' ' << cell->level() << ' '
4452  << cell->material_id() << '\n';
4453  }
4454  out << '\n'; // end of third line
4455  }
4456  else if (cell->reference_cell() == ReferenceCells::Pyramid)
4457  {
4458  // Draw the pyramid as a collections of two lines.
4459  for (const unsigned int v : {0, 1, 3, 2, 0, 4, 1})
4460  {
4461  out << cell->vertex(v) << ' ' << cell->level() << ' '
4462  << cell->material_id() << '\n';
4463  }
4464  out << '\n'; // end of first line
4465 
4466  for (const unsigned int v : {2, 4, 3})
4467  {
4468  out << cell->vertex(v) << ' ' << cell->level() << ' '
4469  << cell->material_id() << '\n';
4470  }
4471  out << '\n'; // end of second line
4472  }
4473  else
4474  Assert(false, ExcNotImplemented());
4475  }
4476  else // need to handle curved boundaries
4477  {
4478  Assert(cell->reference_cell() == ReferenceCells::Hexahedron,
4479  ExcNotImplemented());
4480  for (const unsigned int face_no :
4482  {
4483  const typename ::Triangulation<dim,
4484  spacedim>::face_iterator
4485  face = cell->face(face_no);
4486 
4487  if (face->at_boundary() &&
4488  gnuplot_flags.write_additional_boundary_lines)
4489  {
4490  const unsigned int offset = face_no * n_points * n_points;
4491  for (unsigned int i = 0; i < n_points - 1; ++i)
4492  for (unsigned int j = 0; j < n_points - 1; ++j)
4493  {
4494  const Point<spacedim> p0 =
4495  mapping->transform_unit_to_real_cell(
4496  cell,
4497  q_projector->point(offset + i * n_points + j));
4498  out << p0 << ' ' << cell->level() << ' '
4499  << cell->material_id() << '\n';
4500  out << (mapping->transform_unit_to_real_cell(
4501  cell,
4502  q_projector->point(
4503  offset + (i + 1) * n_points + j)))
4504  << ' ' << cell->level() << ' '
4505  << cell->material_id() << '\n';
4506  out << (mapping->transform_unit_to_real_cell(
4507  cell,
4508  q_projector->point(
4509  offset + (i + 1) * n_points + j + 1)))
4510  << ' ' << cell->level() << ' '
4511  << cell->material_id() << '\n';
4512  out << (mapping->transform_unit_to_real_cell(
4513  cell,
4514  q_projector->point(offset + i * n_points +
4515  j + 1)))
4516  << ' ' << cell->level() << ' '
4517  << cell->material_id() << '\n';
4518  // and the first point again
4519  out << p0 << ' ' << cell->level() << ' '
4520  << cell->material_id() << '\n';
4521  out << '\n' << '\n';
4522  }
4523  }
4524  else
4525  {
4526  for (unsigned int l = 0;
4527  l < GeometryInfo<dim>::lines_per_face;
4528  ++l)
4529  {
4530  const typename ::Triangulation<dim, spacedim>::
4531  line_iterator line = face->line(l);
4532 
4533  const Point<spacedim> &v0 = line->vertex(0),
4534  &v1 = line->vertex(1);
4535  if (line->at_boundary() ||
4536  gnuplot_flags.curved_inner_cells)
4537  {
4538  // Save the points on each face to a vector and
4539  // then try to remove colinear points that won't
4540  // show up in the generated plot.
4541  std::vector<Point<spacedim>> line_points;
4542  // transform_real_to_unit_cell could be replaced
4543  // by using QProjector<dim>::project_to_line
4544  // which is not yet implemented
4545  const Point<spacedim>
4546  u0 = mapping->transform_real_to_unit_cell(cell,
4547  v0),
4548  u1 = mapping->transform_real_to_unit_cell(cell,
4549  v1);
4550  for (unsigned int i = 0; i < n_points; ++i)
4551  line_points.push_back(
4552  mapping->transform_unit_to_real_cell(
4553  cell,
4554  (1 - boundary_points[i][0]) * u0 +
4555  boundary_points[i][0] * u1));
4556  internal::remove_colinear_points(line_points);
4557  for (const Point<spacedim> &point : line_points)
4558  out << point << ' ' << cell->level() << ' '
4559  << static_cast<unsigned int>(
4560  cell->material_id())
4561  << '\n';
4562  }
4563  else
4564  out << v0 << ' ' << cell->level() << ' '
4565  << cell->material_id() << '\n'
4566  << v1 << ' ' << cell->level() << ' '
4567  << cell->material_id() << '\n';
4568 
4569  out << '\n' << '\n';
4570  }
4571  }
4572  }
4573  }
4574  }
4575 
4576  // make sure everything now gets to disk
4577  out.flush();
4578 
4579  AssertThrow(out.fail() == false, ExcIO());
4580  }
4581  } // namespace
4582 } // namespace internal
4583 
4584 
4585 
4586 template <int dim, int spacedim>
4587 void
4589  std::ostream & out,
4590  const Mapping<dim, spacedim> * mapping) const
4591 {
4592  internal::write_gnuplot(tria, out, mapping, gnuplot_flags);
4593 }
4594 
4595 
4596 
4597 namespace internal
4598 {
4599  namespace
4600  {
4601  struct LineEntry
4602  {
4605  bool colorize;
4606  unsigned int level;
4607  LineEntry(const Point<2> & f,
4608  const Point<2> & s,
4609  const bool c,
4610  const unsigned int l)
4611  : first(f)
4612  , second(s)
4613  , colorize(c)
4614  , level(l)
4615  {}
4616  };
4617 
4618 
4619  void
4620  write_eps(const ::Triangulation<1> &,
4621  std::ostream &,
4622  const Mapping<1> *,
4623  const GridOutFlags::Eps<2> &,
4624  const GridOutFlags::Eps<3> &)
4625  {
4626  Assert(false, ExcNotImplemented());
4627  }
4628 
4629  void
4630  write_eps(const ::Triangulation<1, 2> &,
4631  std::ostream &,
4632  const Mapping<1, 2> *,
4633  const GridOutFlags::Eps<2> &,
4634  const GridOutFlags::Eps<3> &)
4635  {
4636  Assert(false, ExcNotImplemented());
4637  }
4638 
4639  void
4640  write_eps(const ::Triangulation<1, 3> &,
4641  std::ostream &,
4642  const Mapping<1, 3> *,
4643  const GridOutFlags::Eps<2> &,
4644  const GridOutFlags::Eps<3> &)
4645  {
4646  Assert(false, ExcNotImplemented());
4647  }
4648 
4649  void
4650  write_eps(const ::Triangulation<2, 3> &,
4651  std::ostream &,
4652  const Mapping<2, 3> *,
4653  const GridOutFlags::Eps<2> &,
4654  const GridOutFlags::Eps<3> &)
4655  {
4656  Assert(false, ExcNotImplemented());
4657  }
4658 
4659 
4660 
4661  template <int dim, int spacedim>
4662  void
4663  write_eps(const ::Triangulation<dim, spacedim> &tria,
4664  std::ostream & out,
4665  const Mapping<dim, spacedim> * mapping,
4666  const GridOutFlags::Eps<2> & eps_flags_2,
4667  const GridOutFlags::Eps<3> & eps_flags_3)
4668  {
4669  using LineList = std::list<LineEntry>;
4670 
4671  // We should never get here in 1D since this function is overloaded for
4672  // all dim == 1 cases.
4673  Assert(dim == 2 || dim == 3, ExcInternalError());
4674 
4675  // Copy, with an object slice, something containing the flags common to
4676  // all dimensions in order to avoid the recurring distinctions between
4677  // the different eps_flags present.
4678  const GridOutFlags::EpsFlagsBase eps_flags_base =
4679  dim == 2 ?
4680  static_cast<const GridOutFlags::EpsFlagsBase &>(eps_flags_2) :
4681  static_cast<const GridOutFlags::EpsFlagsBase &>(eps_flags_3);
4682 
4683  AssertThrow(out.fail() == false, ExcIO());
4684  const unsigned int n_points = eps_flags_base.n_boundary_face_points;
4685 
4686  // make up a list of lines by which
4687  // we will construct the triangulation
4688  //
4689  // this part unfortunately is a bit
4690  // dimension dependent, so we have to
4691  // treat every dimension different.
4692  // however, by directly producing
4693  // the lines to be printed, i.e. their
4694  // 2d images, we can later do the
4695  // actual output dimension independent
4696  // again
4697  LineList line_list;
4698 
4699  switch (dim)
4700  {
4701  case 1:
4702  {
4703  Assert(false, ExcInternalError());
4704  break;
4705  }
4706 
4707  case 2:
4708  {
4709  for (const auto &cell : tria.active_cell_iterators())
4710  for (const unsigned int line_no : cell->line_indices())
4711  {
4712  typename ::Triangulation<dim, spacedim>::line_iterator
4713  line = cell->line(line_no);
4714 
4715  // first treat all
4716  // interior lines and
4717  // make up a list of
4718  // them. if curved
4719  // lines shall not be
4720  // supported (i.e. no
4721  // mapping is
4722  // provided), then also
4723  // treat all other
4724  // lines
4725  if (!line->has_children() &&
4726  (mapping == nullptr || !line->at_boundary()))
4727  // one would expect
4728  // make_pair(line->vertex(0),
4729  // line->vertex(1))
4730  // here, but that is
4731  // not dimension
4732  // independent, since
4733  // vertex(i) is
4734  // Point<dim>, but we
4735  // want a Point<2>.
4736  // in fact, whenever
4737  // we're here, the
4738  // vertex is a
4739  // Point<dim>, but
4740  // the compiler does
4741  // not know
4742  // this. hopefully,
4743  // the compiler will
4744  // optimize away this
4745  // little kludge
4746  line_list.emplace_back(
4747  Point<2>(line->vertex(0)(0), line->vertex(0)(1)),
4748  Point<2>(line->vertex(1)(0), line->vertex(1)(1)),
4749  line->user_flag_set(),
4750  cell->level());
4751  }
4752 
4753  // next if we are to treat
4754  // curved boundaries
4755  // specially, then add lines
4756  // to the list consisting of
4757  // pieces of the boundary
4758  // lines
4759  if (mapping != nullptr)
4760  {
4761  // to do so, first
4762  // generate a sequence of
4763  // points on a face and
4764  // project them onto the
4765  // faces of a unit cell
4766  std::vector<Point<dim - 1>> boundary_points(n_points);
4767 
4768  for (unsigned int i = 0; i < n_points; ++i)
4769  boundary_points[i](0) = 1. * (i + 1) / (n_points + 1);
4770 
4771  Quadrature<dim - 1> quadrature(boundary_points);
4772  Quadrature<dim> q_projector(
4774 
4775  // next loop over all
4776  // boundary faces and
4777  // generate the info from
4778  // them
4779  for (const auto &cell : tria.active_cell_iterators())
4780  for (const unsigned int face_no :
4782  {
4783  const typename ::Triangulation<dim, spacedim>::
4784  face_iterator face = cell->face(face_no);
4785 
4786  if (face->at_boundary())
4787  {
4788  Point<dim> p0_dim(face->vertex(0));
4789  Point<2> p0(p0_dim(0), p0_dim(1));
4790 
4791  // loop over
4792  // all pieces
4793  // of the line
4794  // and generate
4795  // line-lets
4796  const unsigned int offset = face_no * n_points;
4797  for (unsigned int i = 0; i < n_points; ++i)
4798  {
4799  const Point<dim> p1_dim(
4800  mapping->transform_unit_to_real_cell(
4801  cell, q_projector.point(offset + i)));
4802  const Point<2> p1(p1_dim(0), p1_dim(1));
4803 
4804  line_list.emplace_back(p0,
4805  p1,
4806  face->user_flag_set(),
4807  cell->level());
4808  p0 = p1;
4809  }
4810 
4811  // generate last piece
4812  const Point<dim> p1_dim(face->vertex(1));
4813  const Point<2> p1(p1_dim(0), p1_dim(1));
4814  line_list.emplace_back(p0,
4815  p1,
4816  face->user_flag_set(),
4817  cell->level());
4818  }
4819  }
4820  }
4821 
4822  break;
4823  }
4824 
4825  case 3:
4826  {
4827  // curved boundary output
4828  // presently not supported
4829  Assert(mapping == nullptr, ExcNotImplemented());
4830 
4831  // loop over all lines and compute their
4832  // projection on the plane perpendicular
4833  // to the direction of sight
4834 
4835  // direction of view equals the unit
4836  // vector of the position of the
4837  // spectator to the origin.
4838  //
4839  // we chose here the viewpoint as in
4840  // gnuplot as default.
4841  //
4842  // TODO:[WB] Fix a potential problem with viewing angles in 3d Eps
4843  // GridOut
4844  // note: the following might be wrong
4845  // if one of the base vectors below
4846  // is in direction of the viewer, but
4847  // I am too tired at present to fix
4848  // this
4849  const double pi = numbers::PI;
4850  const double z_angle = eps_flags_3.azimut_angle;
4851  const double turn_angle = eps_flags_3.turn_angle;
4852  const Point<dim> view_direction(
4853  -std::sin(z_angle * 2. * pi / 360.) *
4854  std::sin(turn_angle * 2. * pi / 360.),
4855  +std::sin(z_angle * 2. * pi / 360.) *
4856  std::cos(turn_angle * 2. * pi / 360.),
4857  -std::cos(z_angle * 2. * pi / 360.));
4858 
4859  // decide about the two unit vectors
4860  // in this plane. we chose the first one
4861  // to be the projection of the z-axis
4862  // to this plane
4863  const Tensor<1, dim> vector1 =
4864  Point<dim>(0, 0, 1) -
4865  ((Point<dim>(0, 0, 1) * view_direction) * view_direction);
4866  const Tensor<1, dim> unit_vector1 = vector1 / vector1.norm();
4867 
4868  // now the third vector is fixed. we
4869  // chose the projection of a more or
4870  // less arbitrary vector to the plane
4871  // perpendicular to the first one
4872  const Tensor<1, dim> vector2 =
4873  (Point<dim>(1, 0, 0) -
4874  ((Point<dim>(1, 0, 0) * view_direction) * view_direction) -
4875  ((Point<dim>(1, 0, 0) * unit_vector1) * unit_vector1));
4876  const Tensor<1, dim> unit_vector2 = vector2 / vector2.norm();
4877 
4878 
4879  for (const auto &cell : tria.active_cell_iterators())
4880  for (const unsigned int line_no : cell->line_indices())
4881  {
4882  typename ::Triangulation<dim, spacedim>::line_iterator
4883  line = cell->line(line_no);
4884  line_list.emplace_back(
4885  Point<2>(line->vertex(0) * unit_vector2,
4886  line->vertex(0) * unit_vector1),
4887  Point<2>(line->vertex(1) * unit_vector2,
4888  line->vertex(1) * unit_vector1),
4889  line->user_flag_set(),
4890  cell->level());
4891  }
4892 
4893  break;
4894  }
4895 
4896  default:
4897  Assert(false, ExcNotImplemented());
4898  }
4899 
4900 
4901 
4902  // find out minimum and maximum x and
4903  // y coordinates to compute offsets
4904  // and scaling factors
4905  double x_min = tria.begin_active()->vertex(0)(0);
4906  double x_max = x_min;
4907  double y_min = tria.begin_active()->vertex(0)(1);
4908  double y_max = y_min;
4909  unsigned int max_level = line_list.begin()->level;
4910 
4911  for (LineList::const_iterator line = line_list.begin();
4912  line != line_list.end();
4913  ++line)
4914  {
4915  x_min = std::min(x_min, line->first(0));
4916  x_min = std::min(x_min, line->second(0));
4917 
4918  x_max = std::max(x_max, line->first(0));
4919  x_max = std::max(x_max, line->second(0));
4920 
4921  y_min = std::min(y_min, line->first(1));
4922  y_min = std::min(y_min, line->second(1));
4923 
4924  y_max = std::max(y_max, line->first(1));
4925  y_max = std::max(y_max, line->second(1));
4926 
4927  max_level = std::max(max_level, line->level);
4928  }
4929 
4930  // scale in x-direction such that
4931  // in the output 0 <= x <= 300.
4932  // don't scale in y-direction to
4933  // preserve the shape of the
4934  // triangulation
4935  const double scale =
4936  (eps_flags_base.size /
4937  (eps_flags_base.size_type == GridOutFlags::EpsFlagsBase::width ?
4938  x_max - x_min :
4939  y_min - y_max));
4940 
4941 
4942  // now write preamble
4943  {
4944  // block this to have local
4945  // variables destroyed after
4946  // use
4947  std::time_t time1 = std::time(nullptr);
4948  std::tm * time = std::localtime(&time1);
4949  out << "%!PS-Adobe-2.0 EPSF-1.2" << '\n'
4950  << "%%Title: deal.II Output" << '\n'
4951  << "%%Creator: the deal.II library" << '\n'
4952  << "%%Creation Date: " << time->tm_year + 1900 << "/"
4953  << time->tm_mon + 1 << "/" << time->tm_mday << " - "
4954  << time->tm_hour << ":" << std::setw(2) << time->tm_min << ":"
4955  << std::setw(2) << time->tm_sec << '\n'
4956  << "%%BoundingBox: "
4957  // lower left corner
4958  << "0 0 "
4959  // upper right corner
4960  << static_cast<unsigned int>(
4961  std::floor(((x_max - x_min) * scale) + 1))
4962  << ' '
4963  << static_cast<unsigned int>(
4964  std::floor(((y_max - y_min) * scale) + 1))
4965  << '\n';
4966 
4967  // define some abbreviations to keep
4968  // the output small:
4969  // m=move turtle to
4970  // x=execute line stroke
4971  // b=black pen
4972  // r=red pen
4973  out << "/m {moveto} bind def" << '\n'
4974  << "/x {lineto stroke} bind def" << '\n'
4975  << "/b {0 0 0 setrgbcolor} def" << '\n'
4976  << "/r {1 0 0 setrgbcolor} def" << '\n';
4977 
4978  // calculate colors for level
4979  // coloring; level 0 is black,
4980  // other levels are blue
4981  // ... red
4982  if (eps_flags_base.color_lines_level)
4983  out << "/l { neg " << (max_level) << " add "
4984  << (0.66666 / std::max(1U, (max_level - 1)))
4985  << " mul 1 0.8 sethsbcolor} def" << '\n';
4986 
4987  // in 2d, we can also plot cell
4988  // and vertex numbers, but this
4989  // requires a somewhat more
4990  // lengthy preamble. please
4991  // don't ask me what most of
4992  // this means, it is reverse
4993  // engineered from what GNUPLOT
4994  // uses in its output
4995  if ((dim == 2) && (eps_flags_2.write_cell_numbers ||
4996  eps_flags_2.write_vertex_numbers))
4997  {
4998  out
4999  << ("/R {rmoveto} bind def\n"
5000  "/Symbol-Oblique /Symbol findfont [1 0 .167 1 0 0] makefont\n"
5001  "dup length dict begin {1 index /FID eq {pop pop} {def} ifelse} forall\n"
5002  "currentdict end definefont\n"
5003  "/MFshow {{dup dup 0 get findfont exch 1 get scalefont setfont\n"
5004  "[ currentpoint ] exch dup 2 get 0 exch rmoveto dup dup 5 get exch 4 get\n"
5005  "{show} {stringwidth pop 0 rmoveto}ifelse dup 3 get\n"
5006  "{2 get neg 0 exch rmoveto pop} {pop aload pop moveto}ifelse} forall} bind def\n"
5007  "/MFwidth {0 exch {dup 3 get{dup dup 0 get findfont exch 1 get scalefont setfont\n"
5008  "5 get stringwidth pop add}\n"
5009  "{pop} ifelse} forall} bind def\n"
5010  "/MCshow { currentpoint stroke m\n"
5011  "exch dup MFwidth -2 div 3 -1 roll R MFshow } def\n")
5012  << '\n';
5013  }
5014 
5015  out << "%%EndProlog" << '\n' << '\n';
5016 
5017  // set fine lines
5018  out << eps_flags_base.line_width << " setlinewidth" << '\n';
5019  }
5020 
5021  // now write the lines
5022  const Point<2> offset(x_min, y_min);
5023 
5024  for (LineList::const_iterator line = line_list.begin();
5025  line != line_list.end();
5026  ++line)
5027  if (eps_flags_base.color_lines_level && (line->level > 0))
5028  out << line->level << " l " << (line->first - offset) * scale << " m "
5029  << (line->second - offset) * scale << " x" << '\n';
5030  else
5031  out << ((line->colorize && eps_flags_base.color_lines_on_user_flag) ?
5032  "r " :
5033  "b ")
5034  << (line->first - offset) * scale << " m "
5035  << (line->second - offset) * scale << " x" << '\n';
5036 
5037  // finally write the cell numbers
5038  // in 2d, if that is desired
5039  if ((dim == 2) && (eps_flags_2.write_cell_numbers == true))
5040  {
5041  out << "(Helvetica) findfont 140 scalefont setfont" << '\n';
5042 
5043  for (const auto &cell : tria.active_cell_iterators())
5044  {
5045  out << (cell->center()(0) - offset(0)) * scale << ' '
5046  << (cell->center()(1) - offset(1)) * scale << " m" << '\n'
5047  << "[ [(Helvetica) 12.0 0.0 true true (";
5048  if (eps_flags_2.write_cell_number_level)
5049  out << cell;
5050  else
5051  out << cell->index();
5052 
5053  out << ")] "
5054  << "] -6 MCshow" << '\n';
5055  }
5056  }
5057 
5058  // and the vertex numbers
5059  if ((dim == 2) && (eps_flags_2.write_vertex_numbers == true))
5060  {
5061  out << "(Helvetica) findfont 140 scalefont setfont" << '\n';
5062 
5063  // have a list of those
5064  // vertices which we have
5065  // already tracked, to avoid
5066  // doing this multiply
5067  std::set<unsigned int> treated_vertices;
5068  for (const auto &cell : tria.active_cell_iterators())
5069  for (const unsigned int vertex_no : cell->vertex_indices())
5070  if (treated_vertices.find(cell->vertex_index(vertex_no)) ==
5071  treated_vertices.end())
5072  {
5073  treated_vertices.insert(cell->vertex_index(vertex_no));
5074 
5075  out << (cell->vertex(vertex_no)(0) - offset(0)) * scale << ' '
5076  << (cell->vertex(vertex_no)(1) - offset(1)) * scale
5077  << " m" << '\n'
5078  << "[ [(Helvetica) 10.0 0.0 true true ("
5079  << cell->vertex_index(vertex_no) << ")] "
5080  << "] -6 MCshow" << '\n';
5081  }
5082  }
5083 
5084  out << "showpage" << '\n';
5085 
5086  // make sure everything now gets to
5087  // disk
5088  out.flush();
5089 
5090  AssertThrow(out.fail() == false, ExcIO());
5091  }
5092  } // namespace
5093 } // namespace internal
5094 
5095 
5096 template <int dim, int spacedim>
5097 void
5099  std::ostream & out,
5100  const Mapping<dim, spacedim> * mapping) const
5101 {
5103 }
5104 
5105 
5106 template <int dim, int spacedim>
5107 void
5109  std::ostream & out,
5110  const OutputFormat output_format,
5111  const Mapping<dim, spacedim> * mapping) const
5112 {
5113  switch (output_format)
5114  {
5115  case none:
5116  return;
5117 
5118  case dx:
5119  write_dx(tria, out);
5120  return;
5121 
5122  case ucd:
5123  write_ucd(tria, out);
5124  return;
5125 
5126  case gnuplot:
5127  write_gnuplot(tria, out, mapping);
5128  return;
5129 
5130  case eps:
5131  write_eps(tria, out, mapping);
5132  return;
5133 
5134  case xfig:
5135  write_xfig(tria, out, mapping);
5136  return;
5137 
5138  case msh:
5139  write_msh(tria, out);
5140  return;
5141 
5142  case svg:
5143  write_svg(tria, out);
5144  return;
5145 
5146  case mathgl:
5147  write_mathgl(tria, out);
5148  return;
5149 
5150  case vtk:
5151  write_vtk(tria, out);
5152  return;
5153 
5154  case vtu:
5155  write_vtu(tria, out);
5156  return;
5157  }
5158 
5159  Assert(false, ExcInternalError());
5160 }
5161 
5162 
5163 template <int dim, int spacedim>
5164 void
5166  std::ostream & out,
5167  const Mapping<dim, spacedim> * mapping) const
5168 {
5169  write(tria, out, default_format, mapping);
5170 }
5171 
5172 
5173 // explicit instantiations
5174 #include "grid_out.inst"
5175 
5176 
void attach_triangulation(const Triangulation< dim, spacedim > &)
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation={})
virtual void build_patches(const unsigned int n_subdivisions=0)
Definition: data_out.cc:1064
unsigned int n_boundary_faces(const Triangulation< dim, spacedim > &tria) const
Definition: grid_out.cc:3775
GridOutFlags::Vtu vtu_flags
Definition: grid_out.h:1615
GridOutFlags::Eps< 2 > eps_flags_2
Definition: grid_out.h:1584
unsigned int write_ucd_lines(const Triangulation< dim, spacedim > &tria, const unsigned int next_element_index, std::ostream &out) const
Definition: grid_out.cc:4096
void parse_parameters(ParameterHandler &param)
Definition: grid_out.cc:700
void write_svg(const Triangulation< 2, 2 > &tria, std::ostream &out) const
Definition: grid_out.cc:1695
void set_flags(const GridOutFlags::DX &flags)
Definition: grid_out.cc:471
GridOutFlags::Eps< 1 > eps_flags_1
Definition: grid_out.h:1578
GridOutFlags::XFig xfig_flags
Definition: grid_out.h:1595
unsigned int write_ucd_faces(const Triangulation< dim, spacedim > &tria, const unsigned int next_element_index, std::ostream &out) const
Definition: grid_out.cc:4054
void write_mathgl(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition: grid_out.cc:2978
std::string default_suffix() const
Definition: grid_out.cc:593
static void declare_parameters(ParameterHandler &param)
Definition: grid_out.cc:651
static OutputFormat parse_output_format(const std::string &format_name)
Definition: grid_out.cc:601
void write_vtk(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition: grid_out.cc:3301
GridOutFlags::Gnuplot gnuplot_flags
Definition: grid_out.h:1572
void write_msh(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition: grid_out.cc:1019
GridOut()
Definition: grid_out.cc:465
unsigned int write_msh_lines(const Triangulation< dim, spacedim > &tria, const unsigned int next_element_index, std::ostream &out) const
Definition: grid_out.cc:3934
void write_eps(const Triangulation< dim, spacedim > &tria, std::ostream &out, const Mapping< dim, spacedim > *mapping=nullptr) const
Definition: grid_out.cc:5098
static std::string get_output_format_names()
Definition: grid_out.cc:644
GridOutFlags::Eps< 3 > eps_flags_3
Definition: grid_out.h:1590
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:5108
void write_vtu(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition: grid_out.cc:3540
GridOutFlags::Ucd ucd_flags
Definition: grid_out.h:1566
void write_dx(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition: grid_out.cc:781
std::size_t memory_consumption() const
Definition: grid_out.cc:746
unsigned int n_boundary_lines(const Triangulation< dim, spacedim > &tria) const
Definition: grid_out.cc:3791
GridOutFlags::DX dx_flags
Definition: grid_out.h:1554
unsigned int write_msh_faces(const Triangulation< dim, spacedim > &tria, const unsigned int next_element_index, std::ostream &out) const
Definition: grid_out.cc:3895
void write_ucd(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition: grid_out.cc:1119
GridOutFlags::Svg svg_flags
Definition: grid_out.h:1600
OutputFormat default_format
Definition: grid_out.h:1549
void write_xfig(const Triangulation< dim, spacedim > &tria, std::ostream &out, const Mapping< dim, spacedim > *mapping=nullptr) const
Definition: grid_out.cc:1236
void write_gnuplot(const Triangulation< dim, spacedim > &tria, std::ostream &out, const Mapping< dim, spacedim > *mapping=nullptr) const
Definition: grid_out.cc:4588
GridOutFlags::Vtk vtk_flags
Definition: grid_out.h:1610
GridOutFlags::Msh msh_flags
Definition: grid_out.h:1560
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:3587
OutputFormat
Definition: grid_out.h:1001
@ vtk
write() calls write_vtk()
Definition: grid_out.h:1021
@ eps
write() calls write_eps()
Definition: grid_out.h:1009
@ msh
write() calls write_msh()
Definition: grid_out.h:1015
@ xfig
write() calls write_xfig()
Definition: grid_out.h:1013
@ dx
write() calls write_dx()
Definition: grid_out.h:1005
@ ucd
write() calls write_ucd()
Definition: grid_out.h:1011
@ gnuplot
write() calls write_gnuplot()
Definition: grid_out.h:1007
@ mathgl
write() calls write_mathgl()
Definition: grid_out.h:1019
@ svg
write() calls write_svg()
Definition: grid_out.h:1017
@ none
Do nothing in write()
Definition: grid_out.h:1003
@ vtu
write() calls write_vtu()
Definition: grid_out.h:1023
GridOutFlags::MathGL mathgl_flags
Definition: grid_out.h:1605
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const =0
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const =0
long int get_integer(const std::string &entry_string) const
bool get_bool(const std::string &entry_name) const
void declare_entry(const std::string &entry, const std::string &default_value, const Patterns::PatternBase &pattern=Patterns::Anything(), const std::string &documentation="", const bool has_to_be_set=false)
std::string get(const std::string &entry_string) const
double get_double(const std::string &entry_name) const
void enter_subsection(const std::string &subsection)
Definition: point.h:111
static Quadrature< dim > project_to_all_faces(const Quadrature< dim - 1 > &quadrature)
Definition: qprojector.h:579
const Point< dim > & point(const unsigned int i) const
unsigned int gmsh_element_type() const
Definition: tensor.h:503
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:2660
numbers::NumberTraits< Number >::real_type norm() const
const std::vector< bool > & get_used_vertices() const
void save_user_flags_line(std::ostream &out) const
void save(Archive &ar, const unsigned int version) const
cell_iterator begin(const unsigned int level=0) const
virtual types::subdomain_id locally_owned_subdomain() const
unsigned int n_active_cells() const
unsigned int n_used_vertices() const
cell_iterator end() const
const std::vector< ReferenceCell > & get_reference_cells() const
const std::vector< Point< spacedim > > & get_vertices() const
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: vector.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:416
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:417
Point< 3 > vertices[4]
Point< 2 > second
Definition: grid_out.cc:4604
bool colorize
Definition: grid_out.cc:4605
Point< 2 > first
Definition: grid_out.cc:4603
unsigned int level
Definition: grid_out.cc:4606
const unsigned int v0
Definition: grid_tools.cc:1000
const unsigned int v1
Definition: grid_tools.cc:1000
IteratorRange< active_cell_iterator > active_cell_iterators() const
IteratorRange< active_face_iterator > active_face_iterators() const
IteratorRange< cell_iterator > cell_iterators() const
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
void write_pvtu_record(std::ostream &out, const std::vector< std::string > &piece_names) const
static ::ExceptionBase & ExcInvalidState()
static ::ExceptionBase & ExcInternalError()
ReferenceCell reference_cell
Table< 2, float > data
static void declare_parameters(ParameterHandler &prm)
#define Assert(cond, exc)
Definition: exceptions.h:1473
Point< spacedim > vertices[GeometryInfo< dim >::vertices_per_cell]
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
void parse_parameters(const ParameterHandler &prm)
unsigned int n_subdivisions
static ::ExceptionBase & ExcIO()
std::array< unsigned int, GeometryInfo< dim >::faces_per_cell > neighbors
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
void write_vtu(const std::vector< Patch< dim, spacedim >> &patches, const std::vector< std::string > &data_names, const std::vector< std::tuple< unsigned int, unsigned int, std::string, DataComponentInterpretation::DataComponentInterpretation >> &nonscalar_data_ranges, const VtkFlags &flags, std::ostream &out)
void write_gnuplot(const std::vector< Patch< dim, spacedim >> &patches, const std::vector< std::string > &data_names, const std::vector< std::tuple< unsigned int, unsigned int, std::string, DataComponentInterpretation::DataComponentInterpretation >> &nonscalar_data_ranges, const GnuplotFlags &flags, std::ostream &out)
void write_vtu_header(std::ostream &out, const VtkFlags &flags)
void write_eps(const std::vector< Patch< 2, spacedim >> &patches, const std::vector< std::string > &data_names, const std::vector< std::tuple< unsigned int, unsigned int, std::string, DataComponentInterpretation::DataComponentInterpretation >> &nonscalar_data_ranges, const EpsFlags &flags, std::ostream &out)
void write_vtu_footer(std::ostream &out)
void write_vtu_main(const std::vector< Patch< dim, spacedim >> &patches, const std::vector< std::string > &data_names, const std::vector< std::tuple< unsigned int, unsigned int, std::string, DataComponentInterpretation::DataComponentInterpretation >> &nonscalar_data_ranges, const VtkFlags &flags, std::ostream &out)
void write_cells(const std::vector< Patch< dim, spacedim >> &patches, StreamType &out)
Expression floor(const Expression &x)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2084
types::global_dof_index size_type
Definition: cuda_kernels.h:45
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:472
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:185
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
constexpr const ReferenceCell Tetrahedron
constexpr const ReferenceCell Quadrilateral
constexpr const ReferenceCell Wedge
constexpr const ReferenceCell Pyramid
constexpr const ReferenceCell Triangle
constexpr const ReferenceCell Hexahedron
constexpr const ReferenceCell Vertex
constexpr const ReferenceCell Line
VectorType::value_type * end(VectorType &V)
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:140
std::string encode_base64(const std::vector< unsigned char > &binary_input)
Definition: utilities.cc:436
std::string compress(const std::string &input)
Definition: utilities.cc:392
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:473
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:13734
const types::boundary_id invalid_boundary_id
Definition: types.h:244
static constexpr double PI
Definition: numbers.h:233
const types::boundary_id internal_face_boundary_id
Definition: types.h:260
const types::subdomain_id artificial_subdomain_id
Definition: types.h:298
static const unsigned int invalid_unsigned_int
Definition: types.h:201
const types::manifold_id flat_manifold_id
Definition: types.h:269
unsigned int manifold_id
Definition: types.h:141
unsigned int material_id
Definition: types.h:152
unsigned int boundary_id
Definition: types.h:129
static void declare_parameters(ParameterHandler &param)
Definition: grid_out.cc:66
bool write_all_faces
Definition: grid_out.h:82
bool write_faces
Definition: grid_out.h:66
bool write_measure
Definition: grid_out.h:76
void parse_parameters(ParameterHandler &param)
Definition: grid_out.cc:94
bool write_diameter
Definition: grid_out.h:71
bool write_cells
Definition: grid_out.h:61
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 void declare_parameters(ParameterHandler &param)
Definition: grid_out.cc:198
unsigned int n_boundary_face_points
Definition: grid_out.h:362
void parse_parameters(ParameterHandler &param)
Definition: grid_out.cc:231
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
void parse_parameters(ParameterHandler &param)
Definition: grid_out.cc:265
void parse_parameters(ParameterHandler &param)
Definition: grid_out.cc:314
void parse_parameters(ParameterHandler &param)
Definition: grid_out.cc:357
unsigned int n_extra_curved_line_points
Definition: grid_out.h:253
void parse_parameters(ParameterHandler &param)
Definition: grid_out.cc:175
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
static void declare_parameters(ParameterHandler &param)
Definition: grid_out.cc:167
bool write_additional_boundary_lines
Definition: grid_out.h:272
void parse_parameters(ParameterHandler &param)
Definition: grid_out.cc:457
static void declare_parameters(ParameterHandler &param)
Definition: grid_out.cc:451
Msh(const bool write_faces=false, const bool write_lines=false)
Definition: grid_out.cc:104
void parse_parameters(ParameterHandler &param)
Definition: grid_out.cc:118
static void declare_parameters(ParameterHandler &param)
Definition: grid_out.cc:110
bool label_material_id
Definition: grid_out.h:769
bool label_level_number
Definition: grid_out.h:759
bool label_boundary_id
Definition: grid_out.h:788
unsigned int height
Definition: grid_out.h:661
bool label_level_subdomain_id
Definition: grid_out.h:779
bool label_subdomain_id
Definition: grid_out.h:774
Background background
Definition: grid_out.h:708
unsigned int line_thickness
Definition: grid_out.h:672
bool convert_level_number_to_height
Definition: grid_out.h:744
float cell_font_scaling
Definition: grid_out.h:755
Coloring coloring
Definition: grid_out.h:740
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
@ level_subdomain_id
Convert the level subdomain id into the cell color.
Definition: grid_out.h:737
@ subdomain_id
Convert the subdomain id into the cell color.
Definition: grid_out.h:735
@ material_id
Convert the material id into the cell color (default)
Definition: grid_out.h:731
@ level_number
Convert the level number into the cell color.
Definition: grid_out.h:733
unsigned int width
Definition: grid_out.h:667
bool label_cell_index
Definition: grid_out.h:764
unsigned int boundary_line_thickness
Definition: grid_out.h:676
float level_height_factor
Definition: grid_out.h:750
Ucd(const bool write_preamble=false, const bool write_faces=false, const bool write_lines=false)
Definition: grid_out.cc:125
static void declare_parameters(ParameterHandler &param)
Definition: grid_out.cc:136
void parse_parameters(ParameterHandler &param)
Definition: grid_out.cc:145
bool output_only_relevant
Definition: grid_out.h:893
bool serialize_triangulation
Definition: grid_out.h:914
unsigned int n_boundary_face_points
Definition: grid_out.h:582
Point< 2 > scaling
Definition: grid_out.h:587
@ level_number
Convert the level into the cell color.
Definition: grid_out.h:564
@ material_id
Convert the material id into the cell color.
Definition: grid_out.h:562
@ level_subdomain_id
Convert the level subdomain id into the cell color.
Definition: grid_out.h:568
@ subdomain_id
Convert the global subdomain id into the cell color.
Definition: grid_out.h:566
void parse_parameters(ParameterHandler &param)
Definition: grid_out.cc:397
static void declare_parameters(ParameterHandler &param)
Definition: grid_out.cc:381
Point< 2 > offset
Definition: grid_out.h:593
enum GridOutFlags::XFig::Coloring color_by
const ::Triangulation< dim, spacedim > & tria