Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-2330-gf6dfc6c370 2025-01-06 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 for (unsigned int i = 0; i < n_procs; ++i)
3694 filenames.push_back(filename_without_extension.substr(pos) +
3695 ".proc" + Utilities::int_to_string(i, 4) +
3696 ".vtu");
3697
3698 const std::string pvtu_filename =
3699 (filename_without_extension + ".pvtu");
3700 std::ofstream pvtu_output(pvtu_filename);
3701
3702 DataOut<dim, spacedim> data_out;
3703 data_out.attach_triangulation(*tr);
3704
3705 // We need a dummy vector with the names of the data values in the
3706 // .vtu files in order that the .pvtu contains reference these values
3707 const Vector<float> dummy_vector(tr->n_active_cells());
3708 data_out.add_data_vector(dummy_vector, "level");
3709 data_out.add_data_vector(dummy_vector, "subdomain");
3710 data_out.add_data_vector(dummy_vector, "level_subdomain");
3711 data_out.add_data_vector(dummy_vector, "proc_writing");
3712
3713 data_out.build_patches();
3714
3715 data_out.write_pvtu_record(pvtu_output, filenames);
3716 }
3717 }
3718
3719 std::ofstream out(new_file);
3720 std::vector<
3721 std::tuple<unsigned int,
3722 unsigned int,
3723 std::string,
3725 vector_data_ranges;
3727 patches, data_names, vector_data_ranges, vtu_flags, out);
3728}
3729
3730
3731
3732unsigned int
3734{
3735 return 0;
3736}
3737
3738unsigned int
3740{
3741 return 0;
3742}
3743
3744
3745unsigned int
3747{
3748 return 0;
3749}
3750
3751unsigned int
3753{
3754 return 0;
3755}
3756
3757unsigned int
3759{
3760 return 0;
3761}
3762
3763unsigned int
3765{
3766 return 0;
3767}
3768
3769unsigned int
3771{
3772 return 0;
3773}
3774
3775unsigned int
3777{
3778 return 0;
3779}
3780
3781
3782
3783template <int dim, int spacedim>
3784unsigned int
3786{
3788 unsigned int n_faces = 0;
3789
3790 for (const auto &face : tria.active_face_iterators())
3791 if ((face->at_boundary()) && (face->boundary_id() != 0))
3792 ++n_faces;
3793
3794 return n_faces;
3795}
3796
3797
3798
3799template <int dim, int spacedim>
3800unsigned int
3802{
3803 // save the user flags for lines so
3804 // we can use these flags to track
3805 // which ones we've already counted
3806 std::vector<bool> line_flags;
3807 const_cast<::Triangulation<dim, spacedim> &>(tria).save_user_flags_line(
3808 line_flags);
3809 const_cast<::Triangulation<dim, spacedim> &>(tria)
3810 .clear_user_flags_line();
3811
3812 unsigned int n_lines = 0;
3813
3814 for (const auto &cell : tria.active_cell_iterators())
3815 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
3816 if (cell->line(l)->at_boundary() && (cell->line(l)->boundary_id() != 0) &&
3817 (cell->line(l)->user_flag_set() == false))
3818 {
3819 ++n_lines;
3820 cell->line(l)->set_user_flag();
3821 }
3822
3823 // at the end, restore the user
3824 // flags for the lines
3825 const_cast<::Triangulation<dim, spacedim> &>(tria).load_user_flags_line(
3826 line_flags);
3827
3828 return n_lines;
3829}
3830
3831
3832
3833unsigned int
3835 const unsigned int next_element_index,
3836 std::ostream &) const
3837{
3838 return next_element_index;
3839}
3840
3841
3842unsigned int
3844 const unsigned int next_element_index,
3845 std::ostream &) const
3846{
3847 return next_element_index;
3848}
3849
3850unsigned int
3852 const unsigned int next_element_index,
3853 std::ostream &) const
3854{
3855 return next_element_index;
3856}
3857
3858
3859unsigned int
3861 const unsigned int next_element_index,
3862 std::ostream &) const
3863{
3864 return next_element_index;
3865}
3866
3867unsigned int
3869 const unsigned int next_element_index,
3870 std::ostream &) const
3871{
3872 return next_element_index;
3873}
3874
3875
3876unsigned int
3878 const unsigned int next_element_index,
3879 std::ostream &) const
3880{
3881 return next_element_index;
3882}
3883
3884
3885unsigned int
3887 const unsigned int next_element_index,
3888 std::ostream &) const
3889{
3890 return next_element_index;
3891}
3892
3893unsigned int
3895 const unsigned int next_element_index,
3896 std::ostream &) const
3897{
3898 return next_element_index;
3899}
3900
3901
3902
3903template <int dim, int spacedim>
3904unsigned int
3906 const unsigned int next_element_index,
3907 std::ostream &out) const
3908{
3909 unsigned int current_element_index = next_element_index;
3910
3911 for (const auto &face : tria.active_face_iterators())
3912 if (face->at_boundary() && (face->boundary_id() != 0))
3913 {
3914 out << current_element_index << ' '
3915 << face->reference_cell().gmsh_element_type() << ' ';
3916 out << static_cast<unsigned int>(face->boundary_id()) << ' '
3917 << static_cast<unsigned int>(face->boundary_id()) << ' '
3918 << face->n_vertices();
3919 // note: vertex numbers are 1-base
3920 for (const unsigned int vertex : face->vertex_indices())
3921 {
3922 if (face->reference_cell() == ReferenceCells::Quadrilateral)
3923 out << ' '
3924 << face->vertex_index(
3926 1;
3927 else if ((face->reference_cell() == ReferenceCells::Triangle) ||
3928 (face->reference_cell() == ReferenceCells::Line))
3929 out << ' ' << face->vertex_index(vertex) + 1;
3930 else
3932 }
3933 out << '\n';
3934
3935 ++current_element_index;
3936 }
3937 return current_element_index;
3938}
3939
3940
3941
3942template <int dim, int spacedim>
3943unsigned int
3945 const unsigned int next_element_index,
3946 std::ostream &out) const
3947{
3948 unsigned int current_element_index = next_element_index;
3949 // save the user flags for lines so
3950 // we can use these flags to track
3951 // which ones we've already taken
3952 // care of
3953 std::vector<bool> line_flags;
3954 const_cast<::Triangulation<dim, spacedim> &>(tria).save_user_flags_line(
3955 line_flags);
3956 const_cast<::Triangulation<dim, spacedim> &>(tria)
3957 .clear_user_flags_line();
3958
3959 for (const auto &cell : tria.active_cell_iterators())
3960 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
3961 if (cell->line(l)->at_boundary() && (cell->line(l)->boundary_id() != 0) &&
3962 (cell->line(l)->user_flag_set() == false))
3963 {
3964 out << next_element_index << ' '
3966 out << static_cast<unsigned int>(cell->line(l)->boundary_id()) << ' '
3967 << static_cast<unsigned int>(cell->line(l)->boundary_id())
3968 << " 2 "; // two vertex indices to follow
3969 // note: vertex numbers are 1-base
3970 for (unsigned int vertex = 0; vertex < 2; ++vertex)
3971 out << ' '
3972 << cell->line(l)->vertex_index(
3974 1;
3975 out << '\n';
3976
3977 // move on to the next line
3978 // but mark the current one
3979 // as taken care of
3980 ++current_element_index;
3981 cell->line(l)->set_user_flag();
3982 }
3983
3984 // at the end, restore the user
3985 // flags for the lines
3986 const_cast<::Triangulation<dim, spacedim> &>(tria).load_user_flags_line(
3987 line_flags);
3988
3989 return current_element_index;
3990}
3991
3992
3993
3994unsigned int
3996 const unsigned int next_element_index,
3997 std::ostream &) const
3998{
3999 return next_element_index;
4000}
4001
4002unsigned int
4004 const unsigned int next_element_index,
4005 std::ostream &) const
4006{
4007 return next_element_index;
4008}
4009
4010unsigned int
4012 const unsigned int next_element_index,
4013 std::ostream &) const
4014{
4015 return next_element_index;
4016}
4017
4018unsigned int
4020 const unsigned int next_element_index,
4021 std::ostream &) const
4022{
4023 return next_element_index;
4024}
4025
4026unsigned int
4028 const unsigned int next_element_index,
4029 std::ostream &) const
4030{
4031 return next_element_index;
4032}
4033
4034
4035unsigned int
4037 const unsigned int next_element_index,
4038 std::ostream &) const
4039{
4040 return next_element_index;
4041}
4042
4043
4044unsigned int
4046 const unsigned int next_element_index,
4047 std::ostream &) const
4048{
4049 return next_element_index;
4050}
4051
4052unsigned int
4054 const unsigned int next_element_index,
4055 std::ostream &) const
4056{
4057 return next_element_index;
4058}
4059
4060
4061
4062template <int dim, int spacedim>
4063unsigned int
4065 const unsigned int next_element_index,
4066 std::ostream &out) const
4067{
4068 unsigned int current_element_index = next_element_index;
4070
4071 for (const auto &face : tria.active_face_iterators())
4072 if (face->at_boundary() && (face->boundary_id() != 0))
4073 {
4074 out << current_element_index << " "
4075 << static_cast<unsigned int>(face->boundary_id()) << " ";
4076 switch (dim)
4077 {
4078 case 2:
4079 out << "line ";
4080 break;
4081 case 3:
4082 out << "quad ";
4083 break;
4084 default:
4086 }
4087 // note: vertex numbers are 1-base
4088 for (unsigned int vertex = 0;
4089 vertex < GeometryInfo<dim>::vertices_per_face;
4090 ++vertex)
4091 out << face->vertex_index(
4093 1
4094 << ' ';
4095 out << '\n';
4096
4097 ++current_element_index;
4098 }
4099 return current_element_index;
4100}
4101
4102
4103
4104template <int dim, int spacedim>
4105unsigned int
4107 const unsigned int next_element_index,
4108 std::ostream &out) const
4109{
4110 unsigned int current_element_index = next_element_index;
4111 // save the user flags for lines so
4112 // we can use these flags to track
4113 // which ones we've already taken
4114 // care of
4115 std::vector<bool> line_flags;
4116 const_cast<::Triangulation<dim, spacedim> &>(tria).save_user_flags_line(
4117 line_flags);
4118 const_cast<::Triangulation<dim, spacedim> &>(tria)
4119 .clear_user_flags_line();
4120
4121 for (const auto &cell : tria.active_cell_iterators())
4122 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
4123 if (cell->line(l)->at_boundary() && (cell->line(l)->boundary_id() != 0) &&
4124 (cell->line(l)->user_flag_set() == false))
4125 {
4126 out << current_element_index << " "
4127 << static_cast<unsigned int>(cell->line(l)->boundary_id())
4128 << " line ";
4129 // note: vertex numbers in ucd format are 1-base
4130 for (unsigned int vertex = 0; vertex < 2; ++vertex)
4131 out << cell->line(l)->vertex_index(
4133 1
4134 << ' ';
4135 out << '\n';
4136
4137 // move on to the next line
4138 // but mark the current one
4139 // as taken care of
4140 ++current_element_index;
4141 cell->line(l)->set_user_flag();
4142 }
4143
4144 // at the end, restore the user
4145 // flags for the lines
4146 const_cast<::Triangulation<dim, spacedim> &>(tria).load_user_flags_line(
4147 line_flags);
4148 return current_element_index;
4149}
4150
4151
4152
4153namespace internal
4154{
4155 namespace
4156 {
4165 template <int spacedim>
4166 void
4167 remove_colinear_points(std::vector<Point<spacedim>> &points)
4168 {
4169 while (points.size() > 2)
4170 {
4171 Tensor<1, spacedim> first_difference = points[1] - points[0];
4172 first_difference /= first_difference.norm();
4173 Tensor<1, spacedim> second_difference = points[2] - points[1];
4174 second_difference /= second_difference.norm();
4175 // If the three points are colinear then remove the middle one.
4176 if ((first_difference - second_difference).norm() < 1e-10)
4177 points.erase(points.begin() + 1);
4178 else
4179 break;
4180 }
4181 }
4182
4183
4184
4185 template <int spacedim>
4186 void
4187 write_gnuplot(const ::Triangulation<1, spacedim> &tria,
4188 std::ostream &out,
4189 const Mapping<1, spacedim> *,
4190 const GridOutFlags::Gnuplot &gnuplot_flags)
4191 {
4192 AssertThrow(out.fail() == false, ExcIO());
4193
4194 for (const auto &cell : tria.active_cell_iterators())
4195 {
4196 if (gnuplot_flags.write_cell_numbers)
4197 out << "# cell " << cell << '\n';
4198
4199 out << cell->vertex(0) << ' ' << cell->level() << ' '
4200 << cell->material_id() << '\n'
4201 << cell->vertex(1) << ' ' << cell->level() << ' '
4202 << cell->material_id() << '\n'
4203 << "\n\n";
4204 }
4205
4206 // make sure everything now gets to
4207 // disk
4208 out.flush();
4209
4210 AssertThrow(out.fail() == false, ExcIO());
4211 }
4212
4213
4214
4215 template <int spacedim>
4216 void
4217 write_gnuplot(const ::Triangulation<2, spacedim> &tria,
4218 std::ostream &out,
4219 const Mapping<2, spacedim> *mapping,
4220 const GridOutFlags::Gnuplot &gnuplot_flags)
4221 {
4222 AssertThrow(out.fail() == false, ExcIO());
4223
4224 const int dim = 2;
4225
4226 const unsigned int n_additional_points =
4227 gnuplot_flags.n_extra_curved_line_points;
4228 const unsigned int n_points = 2 + n_additional_points;
4229
4230 // If we need to plot curved lines then generate a quadrature formula to
4231 // place points via the mapping
4232 Quadrature<dim> q_projector;
4233 std::vector<Point<dim - 1>> boundary_points;
4234 if (mapping != nullptr)
4235 {
4236 boundary_points.resize(n_points);
4237 boundary_points[0][0] = 0;
4238 boundary_points[n_points - 1][0] = 1;
4239 for (unsigned int i = 1; i < n_points - 1; ++i)
4240 boundary_points[i][0] = 1. * i / (n_points - 1);
4241
4242 const std::vector<double> dummy_weights(n_points, 1. / n_points);
4243 const Quadrature<dim - 1> quadrature(boundary_points, dummy_weights);
4244
4246 ReferenceCells::get_hypercube<dim>(), quadrature);
4247 }
4248
4249 static constexpr std::array<unsigned int, 8> local_vertex_numbering = {
4250 {0, 1, 5, 4, 2, 3, 7, 6}};
4251 for (const auto &cell : tria.active_cell_iterators())
4252 {
4253 if (gnuplot_flags.write_cell_numbers)
4254 out << "# cell " << cell << '\n';
4255
4256 if (mapping == nullptr ||
4257 (dim == spacedim ?
4258 (!cell->at_boundary() && !gnuplot_flags.curved_inner_cells) :
4259 // ignore checking for boundary or interior cells in the codim
4260 // 1 case: 'or false' is a no-op
4261 false))
4262 {
4263 // write out the four sides of this cell by putting the four
4264 // points (+ the initial point again) in a row and lifting the
4265 // drawing pencil at the end
4266 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
4267 out << cell->vertex(dim == 3 ?
4268 local_vertex_numbering[i] :
4269 GeometryInfo<dim>::ucd_to_deal[i])
4270 << ' ' << cell->level() << ' ' << cell->material_id()
4271 << '\n';
4272 out << cell->vertex(0) << ' ' << cell->level() << ' '
4273 << cell->material_id() << '\n'
4274 << '\n' // double new line for gnuplot 3d plots
4275 << '\n';
4276 }
4277 else
4278 // cell is at boundary and we are to treat curved boundaries. so
4279 // loop over all faces and draw them as small pieces of lines
4280 {
4281 for (const unsigned int face_no :
4282 GeometryInfo<dim>::face_indices())
4283 {
4284 const typename ::Triangulation<dim,
4285 spacedim>::face_iterator
4286 face = cell->face(face_no);
4287 if (dim != spacedim || face->at_boundary() ||
4288 gnuplot_flags.curved_inner_cells)
4289 {
4290 // Save the points on each face to a vector and then try
4291 // to remove colinear points that won't show up in the
4292 // generated plot.
4293 std::vector<Point<spacedim>> line_points;
4294 // compute offset of quadrature points within set of
4295 // projected points
4296 const auto offset =
4298 cell->reference_cell(),
4299 face_no,
4300 cell->combined_face_orientation(face_no),
4301 n_points);
4302 for (unsigned int i = 0; i < n_points; ++i)
4303 line_points.push_back(
4305 cell, q_projector.point(offset + i)));
4306 internal::remove_colinear_points(line_points);
4307
4308 for (const Point<spacedim> &point : line_points)
4309 out << point << ' ' << cell->level() << ' '
4310 << cell->material_id() << '\n';
4311
4312 out << '\n' << '\n';
4313 }
4314 else
4315 {
4316 // if, however, the face is not at the boundary and we
4317 // don't want to curve anything, then draw it as usual
4318 out << face->vertex(0) << ' ' << cell->level() << ' '
4319 << cell->material_id() << '\n'
4320 << face->vertex(1) << ' ' << cell->level() << ' '
4321 << cell->material_id() << '\n'
4322 << '\n'
4323 << '\n';
4324 }
4325 }
4326 }
4327 }
4328
4329 // make sure everything now gets to disk
4330 out.flush();
4331
4332 AssertThrow(out.fail() == false, ExcIO());
4333 }
4334
4335
4336
4337 template <int spacedim>
4338 void
4339 write_gnuplot(const ::Triangulation<3, spacedim> &tria,
4340 std::ostream &out,
4341 const Mapping<3, spacedim> *mapping,
4342 const GridOutFlags::Gnuplot &gnuplot_flags)
4343 {
4344 AssertThrow(out.fail() == false, ExcIO());
4345
4346 const int dim = 3;
4347
4348 const unsigned int n_additional_points =
4349 gnuplot_flags.n_extra_curved_line_points;
4350 const unsigned int n_points = 2 + n_additional_points;
4351
4352 // If we need to plot curved lines then generate a quadrature formula to
4353 // place points via the mapping
4354 std::unique_ptr<Quadrature<dim>> q_projector;
4355 std::vector<Point<1>> boundary_points;
4356 if (mapping != nullptr)
4357 {
4358 boundary_points.resize(n_points);
4359 boundary_points[0][0] = 0;
4360 boundary_points[n_points - 1][0] = 1;
4361 for (unsigned int i = 1; i < n_points - 1; ++i)
4362 boundary_points[i][0] = 1. * i / (n_points - 1);
4363
4364 const std::vector<double> dummy_weights(n_points, 1. / n_points);
4365 const Quadrature<1> quadrature1d(boundary_points, dummy_weights);
4366
4367 // tensor product of points, only one copy
4368 const QIterated<dim - 1> quadrature(quadrature1d, 1);
4369 q_projector = std::make_unique<Quadrature<dim>>(
4371 ReferenceCells::get_hypercube<dim>(), quadrature));
4372 }
4373
4374 for (const auto &cell : tria.active_cell_iterators())
4375 {
4376 if (gnuplot_flags.write_cell_numbers)
4377 out << "# cell " << cell << '\n';
4378
4379 if (mapping == nullptr || n_points == 2 ||
4380 (!cell->has_boundary_lines() &&
4381 !gnuplot_flags.curved_inner_cells))
4382 {
4383 if (cell->reference_cell() == ReferenceCells::Hexahedron)
4384 {
4385 // front face
4386 out << cell->vertex(0) << ' ' << cell->level() << ' '
4387 << cell->material_id() << '\n'
4388 << cell->vertex(1) << ' ' << cell->level() << ' '
4389 << cell->material_id() << '\n'
4390 << cell->vertex(5) << ' ' << cell->level() << ' '
4391 << cell->material_id() << '\n'
4392 << cell->vertex(4) << ' ' << cell->level() << ' '
4393 << cell->material_id() << '\n'
4394 << cell->vertex(0) << ' ' << cell->level() << ' '
4395 << cell->material_id() << '\n'
4396 << '\n';
4397 // back face
4398 out << cell->vertex(2) << ' ' << cell->level() << ' '
4399 << cell->material_id() << '\n'
4400 << cell->vertex(3) << ' ' << cell->level() << ' '
4401 << cell->material_id() << '\n'
4402 << cell->vertex(7) << ' ' << cell->level() << ' '
4403 << cell->material_id() << '\n'
4404 << cell->vertex(6) << ' ' << cell->level() << ' '
4405 << cell->material_id() << '\n'
4406 << cell->vertex(2) << ' ' << cell->level() << ' '
4407 << cell->material_id() << '\n'
4408 << '\n';
4409
4410 // now for the four connecting lines
4411 out << cell->vertex(0) << ' ' << cell->level() << ' '
4412 << cell->material_id() << '\n'
4413 << cell->vertex(2) << ' ' << cell->level() << ' '
4414 << cell->material_id() << '\n'
4415 << '\n';
4416 out << cell->vertex(1) << ' ' << cell->level() << ' '
4417 << cell->material_id() << '\n'
4418 << cell->vertex(3) << ' ' << cell->level() << ' '
4419 << cell->material_id() << '\n'
4420 << '\n';
4421 out << cell->vertex(5) << ' ' << cell->level() << ' '
4422 << cell->material_id() << '\n'
4423 << cell->vertex(7) << ' ' << cell->level() << ' '
4424 << cell->material_id() << '\n'
4425 << '\n';
4426 out << cell->vertex(4) << ' ' << cell->level() << ' '
4427 << cell->material_id() << '\n'
4428 << cell->vertex(6) << ' ' << cell->level() << ' '
4429 << cell->material_id() << '\n'
4430 << '\n';
4431 }
4432 else if (cell->reference_cell() == ReferenceCells::Tetrahedron)
4433 {
4434 // Draw the tetrahedron as a two collections of lines.
4435 for (const unsigned int v : {0, 1, 2, 0, 3, 2})
4436 {
4437 out << cell->vertex(v) << ' ' << cell->level() << ' '
4438 << cell->material_id() << '\n';
4439 }
4440 out << '\n'; // end of first line
4441
4442 for (const unsigned int v : {3, 1})
4443 {
4444 out << cell->vertex(v) << ' ' << cell->level() << ' '
4445 << cell->material_id() << '\n';
4446 }
4447 out << '\n'; // end of second line
4448 }
4449 else if (cell->reference_cell() == ReferenceCells::Wedge)
4450 {
4451 // Draw the wedge as a collection of three
4452 // lines. The first one wraps around the base,
4453 // goes up to the top, and wraps around that. The
4454 // second and third are just individual lines
4455 // going from base to top.
4456 for (const unsigned int v : {0, 1, 2, 0, 3, 4, 5, 3})
4457 {
4458 out << cell->vertex(v) << ' ' << cell->level() << ' '
4459 << cell->material_id() << '\n';
4460 }
4461 out << '\n'; // end of first line
4462
4463 for (const unsigned int v : {1, 4})
4464 {
4465 out << cell->vertex(v) << ' ' << cell->level() << ' '
4466 << cell->material_id() << '\n';
4467 }
4468 out << '\n'; // end of second line
4469
4470 for (const unsigned int v : {2, 5})
4471 {
4472 out << cell->vertex(v) << ' ' << cell->level() << ' '
4473 << cell->material_id() << '\n';
4474 }
4475 out << '\n'; // end of third line
4476 }
4477 else if (cell->reference_cell() == ReferenceCells::Pyramid)
4478 {
4479 // Draw the pyramid as a collections of two lines.
4480 for (const unsigned int v : {0, 1, 3, 2, 0, 4, 1})
4481 {
4482 out << cell->vertex(v) << ' ' << cell->level() << ' '
4483 << cell->material_id() << '\n';
4484 }
4485 out << '\n'; // end of first line
4486
4487 for (const unsigned int v : {2, 4, 3})
4488 {
4489 out << cell->vertex(v) << ' ' << cell->level() << ' '
4490 << cell->material_id() << '\n';
4491 }
4492 out << '\n'; // end of second line
4493 }
4494 else
4496 }
4497 else // need to handle curved boundaries
4498 {
4499 Assert(cell->reference_cell() == ReferenceCells::Hexahedron,
4501 for (const unsigned int face_no :
4502 GeometryInfo<dim>::face_indices())
4503 {
4504 const typename ::Triangulation<dim,
4505 spacedim>::face_iterator
4506 face = cell->face(face_no);
4507
4508 if (face->at_boundary() &&
4509 gnuplot_flags.write_additional_boundary_lines)
4510 {
4511 const auto offset =
4513 cell->reference_cell(),
4514 face_no,
4515 cell->combined_face_orientation(face_no),
4516 n_points * n_points);
4517 for (unsigned int i = 0; i < n_points - 1; ++i)
4518 for (unsigned int j = 0; j < n_points - 1; ++j)
4519 {
4520 const Point<spacedim> p0 =
4522 cell,
4523 q_projector->point(offset + i * n_points + j));
4524 out << p0 << ' ' << cell->level() << ' '
4525 << cell->material_id() << '\n';
4526 out << (mapping->transform_unit_to_real_cell(
4527 cell,
4528 q_projector->point(
4529 offset + (i + 1) * n_points + j)))
4530 << ' ' << cell->level() << ' '
4531 << cell->material_id() << '\n';
4532 out << (mapping->transform_unit_to_real_cell(
4533 cell,
4534 q_projector->point(
4535 offset + (i + 1) * n_points + j + 1)))
4536 << ' ' << cell->level() << ' '
4537 << cell->material_id() << '\n';
4538 out << (mapping->transform_unit_to_real_cell(
4539 cell,
4540 q_projector->point(offset + i * n_points +
4541 j + 1)))
4542 << ' ' << cell->level() << ' '
4543 << cell->material_id() << '\n';
4544 // and the first point again
4545 out << p0 << ' ' << cell->level() << ' '
4546 << cell->material_id() << '\n';
4547 out << '\n' << '\n';
4548 }
4549 }
4550 else
4551 {
4552 for (unsigned int l = 0;
4553 l < GeometryInfo<dim>::lines_per_face;
4554 ++l)
4555 {
4556 const typename ::Triangulation<dim, spacedim>::
4557 line_iterator line = face->line(l);
4558
4559 const Point<spacedim> &v0 = line->vertex(0),
4560 &v1 = line->vertex(1);
4561 if (line->at_boundary() ||
4562 gnuplot_flags.curved_inner_cells)
4563 {
4564 // Save the points on each face to a vector and
4565 // then try to remove colinear points that won't
4566 // show up in the generated plot.
4567 std::vector<Point<spacedim>> line_points;
4568 // transform_real_to_unit_cell could be replaced
4569 // by using QProjector<dim>::project_to_line
4570 // which is not yet implemented
4571 const Point<spacedim>
4572 u0 = mapping->transform_real_to_unit_cell(cell,
4573 v0),
4574 u1 = mapping->transform_real_to_unit_cell(cell,
4575 v1);
4576 for (unsigned int i = 0; i < n_points; ++i)
4577 line_points.push_back(
4579 cell,
4580 (1 - boundary_points[i][0]) * u0 +
4581 boundary_points[i][0] * u1));
4582 internal::remove_colinear_points(line_points);
4583 for (const Point<spacedim> &point : line_points)
4584 out << point << ' ' << cell->level() << ' '
4585 << static_cast<unsigned int>(
4586 cell->material_id())
4587 << '\n';
4588 }
4589 else
4590 out << v0 << ' ' << cell->level() << ' '
4591 << cell->material_id() << '\n'
4592 << v1 << ' ' << cell->level() << ' '
4593 << cell->material_id() << '\n';
4594
4595 out << '\n' << '\n';
4596 }
4597 }
4598 }
4599 }
4600 }
4601
4602 // make sure everything now gets to disk
4603 out.flush();
4604
4605 AssertThrow(out.fail() == false, ExcIO());
4606 }
4607 } // namespace
4608} // namespace internal
4609
4610
4611
4612template <int dim, int spacedim>
4613void
4615 std::ostream &out,
4616 const Mapping<dim, spacedim> *mapping) const
4617{
4618 internal::write_gnuplot(tria, out, mapping, gnuplot_flags);
4619}
4620
4621
4622
4623namespace internal
4624{
4625 namespace
4626 {
4627 struct LineEntry
4628 {
4632 unsigned int level;
4633 LineEntry(const Point<2> &f,
4634 const Point<2> &s,
4635 const bool c,
4636 const unsigned int l)
4637 : first(f)
4638 , second(s)
4639 , colorize(c)
4640 , level(l)
4641 {}
4642 };
4643
4644
4645 void
4646 write_eps(const ::Triangulation<1> &,
4647 std::ostream &,
4648 const Mapping<1> *,
4649 const GridOutFlags::Eps<2> &,
4650 const GridOutFlags::Eps<3> &)
4651 {
4653 }
4654
4655 void
4656 write_eps(const ::Triangulation<1, 2> &,
4657 std::ostream &,
4658 const Mapping<1, 2> *,
4659 const GridOutFlags::Eps<2> &,
4660 const GridOutFlags::Eps<3> &)
4661 {
4663 }
4664
4665 void
4666 write_eps(const ::Triangulation<1, 3> &,
4667 std::ostream &,
4668 const Mapping<1, 3> *,
4669 const GridOutFlags::Eps<2> &,
4670 const GridOutFlags::Eps<3> &)
4671 {
4673 }
4674
4675 void
4676 write_eps(const ::Triangulation<2, 3> &,
4677 std::ostream &,
4678 const Mapping<2, 3> *,
4679 const GridOutFlags::Eps<2> &,
4680 const GridOutFlags::Eps<3> &)
4681 {
4683 }
4684
4685
4686
4687 template <int dim, int spacedim>
4688 void
4689 write_eps(const ::Triangulation<dim, spacedim> &tria,
4690 std::ostream &out,
4691 const Mapping<dim, spacedim> *mapping,
4692 const GridOutFlags::Eps<2> &eps_flags_2,
4693 const GridOutFlags::Eps<3> &eps_flags_3)
4694 {
4695 using LineList = std::list<LineEntry>;
4696
4697 // We should never get here in 1d since this function is overloaded for
4698 // all dim == 1 cases.
4699 Assert(dim == 2 || dim == 3, ExcInternalError());
4700
4701 // Copy, with an object slice, something containing the flags common to
4702 // all dimensions in order to avoid the recurring distinctions between
4703 // the different eps_flags present.
4704 const GridOutFlags::EpsFlagsBase eps_flags_base =
4705 dim == 2 ?
4706 static_cast<const GridOutFlags::EpsFlagsBase &>(eps_flags_2) :
4707 static_cast<const GridOutFlags::EpsFlagsBase &>(eps_flags_3);
4708
4709 AssertThrow(out.fail() == false, ExcIO());
4710 const unsigned int n_points = eps_flags_base.n_boundary_face_points;
4711
4712 // make up a list of lines by which
4713 // we will construct the triangulation
4714 //
4715 // this part unfortunately is a bit
4716 // dimension dependent, so we have to
4717 // treat every dimension different.
4718 // however, by directly producing
4719 // the lines to be printed, i.e. their
4720 // 2d images, we can later do the
4721 // actual output dimension independent
4722 // again
4723 LineList line_list;
4724
4725 switch (dim)
4726 {
4727 case 1:
4728 {
4730 break;
4731 }
4732
4733 case 2:
4734 {
4735 for (const auto &cell : tria.active_cell_iterators())
4736 for (const unsigned int line_no : cell->line_indices())
4737 {
4738 typename ::Triangulation<dim, spacedim>::line_iterator
4739 line = cell->line(line_no);
4740
4741 // first treat all
4742 // interior lines and
4743 // make up a list of
4744 // them. if curved
4745 // lines shall not be
4746 // supported (i.e. no
4747 // mapping is
4748 // provided), then also
4749 // treat all other
4750 // lines
4751 if (!line->has_children() &&
4752 (mapping == nullptr || !line->at_boundary()))
4753 // one would expect
4754 // make_pair(line->vertex(0),
4755 // line->vertex(1))
4756 // here, but that is
4757 // not dimension
4758 // independent, since
4759 // vertex(i) is
4760 // Point<dim>, but we
4761 // want a Point<2>.
4762 // in fact, whenever
4763 // we're here, the
4764 // vertex is a
4765 // Point<dim>, but
4766 // the compiler does
4767 // not know
4768 // this. hopefully,
4769 // the compiler will
4770 // optimize away this
4771 // little kludge
4772 line_list.emplace_back(
4773 Point<2>(line->vertex(0)[0], line->vertex(0)[1]),
4774 Point<2>(line->vertex(1)[0], line->vertex(1)[1]),
4775 line->user_flag_set(),
4776 cell->level());
4777 }
4778
4779 // next if we are to treat
4780 // curved boundaries
4781 // specially, then add lines
4782 // to the list consisting of
4783 // pieces of the boundary
4784 // lines
4785 if (mapping != nullptr)
4786 {
4787 // to do so, first
4788 // generate a sequence of
4789 // points on a face and
4790 // project them onto the
4791 // faces of a unit cell
4792 std::vector<Point<dim - 1>> boundary_points(n_points);
4793
4794 for (unsigned int i = 0; i < n_points; ++i)
4795 boundary_points[i][0] = 1. * (i + 1) / (n_points + 1);
4796
4797 const Quadrature<dim - 1> quadrature(boundary_points);
4798 const Quadrature<dim> q_projector(
4800 ReferenceCells::get_hypercube<dim>(), quadrature));
4801
4802 // next loop over all
4803 // boundary faces and
4804 // generate the info from
4805 // them
4806 for (const auto &cell : tria.active_cell_iterators())
4807 for (const unsigned int face_no :
4808 GeometryInfo<dim>::face_indices())
4809 {
4810 const typename ::Triangulation<dim, spacedim>::
4811 face_iterator face = cell->face(face_no);
4812
4813 if (face->at_boundary())
4814 {
4815 Point<dim> p0_dim(face->vertex(0));
4816 Point<2> p0(p0_dim[0], p0_dim[1]);
4817
4818 // loop over
4819 // all pieces
4820 // of the line
4821 // and generate
4822 // line-lets
4823 const auto offset =
4825 cell->reference_cell(),
4826 face_no,
4827 cell->combined_face_orientation(face_no),
4828 n_points);
4829 for (unsigned int i = 0; i < n_points; ++i)
4830 {
4831 const Point<dim> p1_dim(
4833 cell, q_projector.point(offset + i)));
4834 const Point<2> p1(p1_dim[0], p1_dim[1]);
4835
4836 line_list.emplace_back(p0,
4837 p1,
4838 face->user_flag_set(),
4839 cell->level());
4840 p0 = p1;
4841 }
4842
4843 // generate last piece
4844 const Point<dim> p1_dim(face->vertex(1));
4845 const Point<2> p1(p1_dim[0], p1_dim[1]);
4846 line_list.emplace_back(p0,
4847 p1,
4848 face->user_flag_set(),
4849 cell->level());
4850 }
4851 }
4852 }
4853
4854 break;
4855 }
4856
4857 case 3:
4858 {
4859 // curved boundary output
4860 // presently not supported
4861 Assert(mapping == nullptr, ExcNotImplemented());
4862
4863 // loop over all lines and compute their
4864 // projection on the plane perpendicular
4865 // to the direction of sight
4866
4867 // direction of view equals the unit
4868 // vector of the position of the
4869 // spectator to the origin.
4870 //
4871 // we chose here the viewpoint as in
4872 // gnuplot as default.
4873 //
4874 // TODO:[WB] Fix a potential problem with viewing angles in 3d Eps
4875 // GridOut
4876 // note: the following might be wrong
4877 // if one of the base vectors below
4878 // is in direction of the viewer, but
4879 // I am too tired at present to fix
4880 // this
4881 const double pi = numbers::PI;
4882 const double z_angle = eps_flags_3.azimut_angle;
4883 const double turn_angle = eps_flags_3.turn_angle;
4884 const Point<dim> view_direction(
4885 -std::sin(z_angle * 2. * pi / 360.) *
4886 std::sin(turn_angle * 2. * pi / 360.),
4887 +std::sin(z_angle * 2. * pi / 360.) *
4888 std::cos(turn_angle * 2. * pi / 360.),
4889 -std::cos(z_angle * 2. * pi / 360.));
4890
4891 // decide about the two unit vectors
4892 // in this plane. we chose the first one
4893 // to be the projection of the z-axis
4894 // to this plane
4895 const Tensor<1, dim> vector1 =
4896 Point<dim>(0, 0, 1) -
4897 ((Point<dim>(0, 0, 1) * view_direction) * view_direction);
4898 const Tensor<1, dim> unit_vector1 = vector1 / vector1.norm();
4899
4900 // now the third vector is fixed. we
4901 // chose the projection of a more or
4902 // less arbitrary vector to the plane
4903 // perpendicular to the first one
4904 const Tensor<1, dim> vector2 =
4905 (Point<dim>(1, 0, 0) -
4906 ((Point<dim>(1, 0, 0) * view_direction) * view_direction) -
4907 ((Point<dim>(1, 0, 0) * unit_vector1) * unit_vector1));
4908 const Tensor<1, dim> unit_vector2 = vector2 / vector2.norm();
4909
4910
4911 for (const auto &cell : tria.active_cell_iterators())
4912 for (const unsigned int line_no : cell->line_indices())
4913 {
4914 typename ::Triangulation<dim, spacedim>::line_iterator
4915 line = cell->line(line_no);
4916 line_list.emplace_back(
4917 Point<2>(line->vertex(0) * unit_vector2,
4918 line->vertex(0) * unit_vector1),
4919 Point<2>(line->vertex(1) * unit_vector2,
4920 line->vertex(1) * unit_vector1),
4921 line->user_flag_set(),
4922 cell->level());
4923 }
4924
4925 break;
4926 }
4927
4928 default:
4930 }
4931
4932
4933
4934 // find out minimum and maximum x and
4935 // y coordinates to compute offsets
4936 // and scaling factors
4937 double x_min = tria.begin_active()->vertex(0)[0];
4938 double x_max = x_min;
4939 double y_min = tria.begin_active()->vertex(0)[1];
4940 double y_max = y_min;
4941 unsigned int max_level = line_list.begin()->level;
4942
4943 for (LineList::const_iterator line = line_list.begin();
4944 line != line_list.end();
4945 ++line)
4946 {
4947 x_min = std::min(x_min, line->first[0]);
4948 x_min = std::min(x_min, line->second[0]);
4949
4950 x_max = std::max(x_max, line->first[0]);
4951 x_max = std::max(x_max, line->second[0]);
4952
4953 y_min = std::min(y_min, line->first[1]);
4954 y_min = std::min(y_min, line->second[1]);
4955
4956 y_max = std::max(y_max, line->first[1]);
4957 y_max = std::max(y_max, line->second[1]);
4958
4959 max_level = std::max(max_level, line->level);
4960 }
4961
4962 // scale in x-direction such that
4963 // in the output 0 <= x <= 300.
4964 // don't scale in y-direction to
4965 // preserve the shape of the
4966 // triangulation
4967 const double scale =
4968 (eps_flags_base.size /
4969 (eps_flags_base.size_type == GridOutFlags::EpsFlagsBase::width ?
4970 x_max - x_min :
4971 y_min - y_max));
4972
4973
4974 // now write preamble
4975 {
4976 // block this to have local
4977 // variables destroyed after
4978 // use
4979 std::time_t time1 = std::time(nullptr);
4980 std::tm *time = std::localtime(&time1);
4981 out << "%!PS-Adobe-2.0 EPSF-1.2" << '\n'
4982 << "%%Title: deal.II Output" << '\n'
4983 << "%%Creator: the deal.II library" << '\n'
4984 << "%%Creation Date: " << time->tm_year + 1900 << "/"
4985 << time->tm_mon + 1 << "/" << time->tm_mday << " - "
4986 << time->tm_hour << ":" << std::setw(2) << time->tm_min << ":"
4987 << std::setw(2) << time->tm_sec << '\n'
4988 << "%%BoundingBox: "
4989 // lower left corner
4990 << "0 0 "
4991 // upper right corner
4992 << static_cast<unsigned int>(
4993 std::floor(((x_max - x_min) * scale) + 1))
4994 << ' '
4995 << static_cast<unsigned int>(
4996 std::floor(((y_max - y_min) * scale) + 1))
4997 << '\n';
4998
4999 // define some abbreviations to keep
5000 // the output small:
5001 // m=move turtle to
5002 // x=execute line stroke
5003 // b=black pen
5004 // r=red pen
5005 out << "/m {moveto} bind def" << '\n'
5006 << "/x {lineto stroke} bind def" << '\n'
5007 << "/b {0 0 0 setrgbcolor} def" << '\n'
5008 << "/r {1 0 0 setrgbcolor} def" << '\n';
5009
5010 // calculate colors for level
5011 // coloring; level 0 is black,
5012 // other levels are blue
5013 // ... red
5014 if (eps_flags_base.color_lines_level)
5015 out << "/l { neg " << (max_level) << " add "
5016 << (0.66666 / std::max(1U, (max_level - 1)))
5017 << " mul 1 0.8 sethsbcolor} def" << '\n';
5018
5019 // in 2d, we can also plot cell
5020 // and vertex numbers, but this
5021 // requires a somewhat more
5022 // lengthy preamble. please
5023 // don't ask me what most of
5024 // this means, it is reverse
5025 // engineered from what GNUPLOT
5026 // uses in its output
5027 if ((dim == 2) && (eps_flags_2.write_cell_numbers ||
5028 eps_flags_2.write_vertex_numbers))
5029 {
5030 out
5031 << ("/R {rmoveto} bind def\n"
5032 "/Symbol-Oblique /Symbol findfont [1 0 .167 1 0 0] makefont\n"
5033 "dup length dict begin {1 index /FID eq {pop pop} {def} ifelse} forall\n"
5034 "currentdict end definefont\n"
5035 "/MFshow {{dup dup 0 get findfont exch 1 get scalefont setfont\n"
5036 "[ currentpoint ] exch dup 2 get 0 exch rmoveto dup dup 5 get exch 4 get\n"
5037 "{show} {stringwidth pop 0 rmoveto}ifelse dup 3 get\n"
5038 "{2 get neg 0 exch rmoveto pop} {pop aload pop moveto}ifelse} forall} bind def\n"
5039 "/MFwidth {0 exch {dup 3 get{dup dup 0 get findfont exch 1 get scalefont setfont\n"
5040 "5 get stringwidth pop add}\n"
5041 "{pop} ifelse} forall} bind def\n"
5042 "/MCshow { currentpoint stroke m\n"
5043 "exch dup MFwidth -2 div 3 -1 roll R MFshow } def\n")
5044 << '\n';
5045 }
5046
5047 out << "%%EndProlog" << '\n' << '\n';
5048
5049 // set fine lines
5050 out << eps_flags_base.line_width << " setlinewidth" << '\n';
5051 }
5052
5053 // now write the lines
5054 const Point<2> offset(x_min, y_min);
5055
5056 for (LineList::const_iterator line = line_list.begin();
5057 line != line_list.end();
5058 ++line)
5059 if (eps_flags_base.color_lines_level && (line->level > 0))
5060 out << line->level << " l " << (line->first - offset) * scale << " m "
5061 << (line->second - offset) * scale << " x" << '\n';
5062 else
5063 out << ((line->colorize && eps_flags_base.color_lines_on_user_flag) ?
5064 "r " :
5065 "b ")
5066 << (line->first - offset) * scale << " m "
5067 << (line->second - offset) * scale << " x" << '\n';
5068
5069 // finally write the cell numbers
5070 // in 2d, if that is desired
5071 if ((dim == 2) && (eps_flags_2.write_cell_numbers == true))
5072 {
5073 out << "(Helvetica) findfont 140 scalefont setfont" << '\n';
5074
5075 for (const auto &cell : tria.active_cell_iterators())
5076 {
5077 out << (cell->center()[0] - offset[0]) * scale << ' '
5078 << (cell->center()[1] - offset[1]) * scale << " m" << '\n'
5079 << "[ [(Helvetica) 12.0 0.0 true true (";
5080 if (eps_flags_2.write_cell_number_level)
5081 out << cell;
5082 else
5083 out << cell->index();
5084
5085 out << ")] "
5086 << "] -6 MCshow" << '\n';
5087 }
5088 }
5089
5090 // and the vertex numbers
5091 if ((dim == 2) && (eps_flags_2.write_vertex_numbers == true))
5092 {
5093 out << "(Helvetica) findfont 140 scalefont setfont" << '\n';
5094
5095 // have a list of those
5096 // vertices which we have
5097 // already tracked, to avoid
5098 // doing this multiply
5099 std::set<unsigned int> treated_vertices;
5100 for (const auto &cell : tria.active_cell_iterators())
5101 for (const unsigned int vertex_no : cell->vertex_indices())
5102 if (treated_vertices.find(cell->vertex_index(vertex_no)) ==
5103 treated_vertices.end())
5104 {
5105 treated_vertices.insert(cell->vertex_index(vertex_no));
5106
5107 out << (cell->vertex(vertex_no)[0] - offset[0]) * scale << ' '
5108 << (cell->vertex(vertex_no)[1] - offset[1]) * scale
5109 << " m" << '\n'
5110 << "[ [(Helvetica) 10.0 0.0 true true ("
5111 << cell->vertex_index(vertex_no) << ")] "
5112 << "] -6 MCshow" << '\n';
5113 }
5114 }
5115
5116 out << "showpage" << '\n';
5117
5118 // make sure everything now gets to
5119 // disk
5120 out.flush();
5121
5122 AssertThrow(out.fail() == false, ExcIO());
5123 }
5124 } // namespace
5125} // namespace internal
5126
5127
5128template <int dim, int spacedim>
5129void
5131 std::ostream &out,
5132 const Mapping<dim, spacedim> *mapping) const
5133{
5134 internal::write_eps(tria, out, mapping, eps_flags_2, eps_flags_3);
5135}
5136
5137
5138template <int dim, int spacedim>
5139void
5141 std::ostream &out,
5142 const OutputFormat output_format,
5143 const Mapping<dim, spacedim> *mapping) const
5144{
5145 switch (output_format)
5146 {
5147 case none:
5148 return;
5149
5150 case dx:
5151 write_dx(tria, out);
5152 return;
5153
5154 case ucd:
5155 write_ucd(tria, out);
5156 return;
5157
5158 case gnuplot:
5159 write_gnuplot(tria, out, mapping);
5160 return;
5161
5162 case eps:
5163 write_eps(tria, out, mapping);
5164 return;
5165
5166 case xfig:
5167 write_xfig(tria, out, mapping);
5168 return;
5169
5170 case msh:
5171 write_msh(tria, out);
5172 return;
5173
5174 case svg:
5175 write_svg(tria, out);
5176 return;
5177
5178 case mathgl:
5179 write_mathgl(tria, out);
5180 return;
5181
5182 case vtk:
5183 write_vtk(tria, out);
5184 return;
5185
5186 case vtu:
5187 write_vtu(tria, out);
5188 return;
5189 }
5190
5192}
5193
5194
5195template <int dim, int spacedim>
5196void
5198 std::ostream &out,
5199 const Mapping<dim, spacedim> *mapping) const
5200{
5201 write(tria, out, default_format, mapping);
5202}
5203
5204
5205// explicit instantiations
5206#include "grid_out.inst"
5207
5208
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:3785
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:4106
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:4064
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:3944
void write_eps(const Triangulation< dim, spacedim > &tria, std::ostream &out, const Mapping< dim, spacedim > *mapping=nullptr) const
Definition grid_out.cc:5130
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:5140
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:3801
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:3905
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:4614
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:318
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:111
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)
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:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
Point< 2 > second
Definition grid_out.cc:4630
bool colorize
Definition grid_out.cc:4631
Point< 2 > first
Definition grid_out.cc:4629
unsigned int level
Definition grid_out.cc:4632
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)
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
std::size_t size
Definition mpi.cc:734
const unsigned int n_procs
Definition mpi.cc:924
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:191
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:92
std::string encode_base64(const std::vector< unsigned char > &binary_input)
Definition utilities.cc:433
std::string compress(const std::string &input)
Definition utilities.cc:389
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:470
const types::boundary_id invalid_boundary_id
Definition types.h:292
static constexpr double PI
Definition numbers.h:254
const types::boundary_id internal_face_boundary_id
Definition types.h:312
const types::subdomain_id artificial_subdomain_id
Definition types.h:362
static const unsigned int invalid_unsigned_int
Definition types.h:220
const types::manifold_id flat_manifold_id
Definition types.h:325
::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:167
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