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