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