deal.II version GIT relicensing-3323-g74e9678af2 2025-05-17 01:50: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\}}\)
Loading...
Searching...
No Matches
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 IteratorType>
3127 void
3128 generate_triangulation_patches(
3129 std::vector<DataOutBase::Patch<dim, spacedim>> &patches,
3130 const IteratorType &begin,
3132 {
3133 // convert each of the active cells into a patch
3134 for (auto cell = begin; cell != end; ++cell)
3135 {
3137 patch.reference_cell = cell->reference_cell();
3138 patch.n_subdivisions = 1;
3139 patch.data.reinit(5, cell->n_vertices());
3140
3141 for (const unsigned int v : cell->vertex_indices())
3142 {
3143 patch.vertices[v] = cell->vertex(v);
3144 patch.data(0, v) = cell->level();
3145 patch.data(1, v) =
3146 static_cast<std::make_signed_t<types::manifold_id>>(
3147 cell->manifold_id());
3148 patch.data(2, v) =
3149 static_cast<std::make_signed_t<types::material_id>>(
3150 cell->material_id());
3151 if (cell->is_active())
3152 patch.data(3, v) =
3153 static_cast<std::make_signed_t<types::subdomain_id>>(
3154 cell->subdomain_id());
3155 else
3156 patch.data(3, v) = -1;
3157 patch.data(4, v) =
3158 static_cast<std::make_signed_t<types::subdomain_id>>(
3159 cell->level_subdomain_id());
3160 }
3161 patches.push_back(patch);
3162 }
3163 }
3164
3165
3166
3167 std::vector<std::string>
3168 triangulation_patch_data_names()
3169 {
3170 std::vector<std::string> v(5);
3171 v[0] = "level";
3172 v[1] = "manifold";
3173 v[2] = "material";
3174 v[3] = "subdomain";
3175 v[4] = "level_subdomain";
3176 return v;
3177 }
3178
3182 std::vector<typename Triangulation<3, 3>::active_line_iterator>
3183 get_boundary_edge_iterators(const Triangulation<3, 3> &tria)
3184 {
3185 std::vector<typename Triangulation<3, 3>::active_line_iterator> res;
3186
3187 std::vector<bool> flags;
3188 tria.save_user_flags_line(flags);
3189 const_cast<Triangulation<3, 3> &>(tria).clear_user_flags_line();
3190
3191 for (auto face : tria.active_face_iterators())
3192 for (const auto l : face->line_indices())
3193 {
3194 const auto line = face->line(l);
3195 if (line->user_flag_set() || line->has_children())
3196 continue;
3197 else
3198 line->set_user_flag();
3199 if (line->at_boundary())
3200 res.emplace_back(line);
3201 }
3202 const_cast<Triangulation<3, 3> &>(tria).load_user_flags_line(flags);
3203 return res;
3204 }
3205
3206
3207
3211 template <int dim, int spacedim>
3212 std::vector<typename Triangulation<dim, spacedim>::active_line_iterator>
3213 get_boundary_edge_iterators(const Triangulation<dim, spacedim> &)
3214 {
3215 return {};
3216 }
3217
3218
3219
3224 std::vector<typename Triangulation<3, 3>::active_line_iterator>
3225 get_relevant_edge_iterators(const Triangulation<3, 3> &tria)
3226 {
3227 std::vector<typename Triangulation<3, 3>::active_line_iterator> res;
3228
3229 std::vector<bool> flags;
3230 tria.save_user_flags_line(flags);
3231 const_cast<Triangulation<3, 3> &>(tria).clear_user_flags_line();
3232
3233 for (auto face : tria.active_face_iterators())
3234 for (const auto l : face->line_indices())
3235 {
3236 const auto line = face->line(l);
3237 if (line->user_flag_set() || line->has_children())
3238 continue;
3239 else
3240 line->set_user_flag();
3241 if (line->manifold_id() != numbers::flat_manifold_id ||
3242 (line->boundary_id() != 0 &&
3243 line->boundary_id() != numbers::invalid_boundary_id))
3244 res.emplace_back(line);
3245 }
3246 const_cast<Triangulation<3, 3> &>(tria).load_user_flags_line(flags);
3247 return res;
3248 }
3249
3250
3254 template <int dim, int spacedim>
3255 std::vector<typename Triangulation<dim, spacedim>::active_line_iterator>
3256 get_relevant_edge_iterators(const Triangulation<dim, spacedim> &)
3257 {
3258 return {};
3259 }
3260
3261
3262
3266 template <int dim, int spacedim>
3267 std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3268 get_boundary_face_iterators(const Triangulation<dim, spacedim> &tria)
3269 {
3270 std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3271 res;
3272 if (dim == 1)
3273 return res;
3274 for (auto face : tria.active_face_iterators())
3275 {
3276 if (face->boundary_id() != numbers::invalid_boundary_id)
3277 res.push_back(face);
3278 }
3279 return res;
3280 }
3281
3282
3283
3288 template <int dim, int spacedim>
3289 std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3290 get_relevant_face_iterators(const Triangulation<dim, spacedim> &tria)
3291 {
3292 std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3293 res;
3294 if (dim == 1)
3295 return res;
3296 for (auto face : tria.active_face_iterators())
3297 {
3298 if (face->manifold_id() != numbers::flat_manifold_id ||
3299 (face->boundary_id() != 0 &&
3300 face->boundary_id() != numbers::invalid_boundary_id))
3301 res.push_back(face);
3302 }
3303 return res;
3304 }
3305} // namespace
3306
3307
3308
3309template <int dim, int spacedim>
3310void
3312 std::ostream &out) const
3313{
3314 AssertThrow(out.fail() == false, ExcIO());
3315
3316 // get the positions of the vertices
3317 const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
3318
3319 const auto n_vertices = vertices.size();
3320
3321 out << "# vtk DataFile Version 3.0\n"
3322 << "Triangulation generated with deal.II\n"
3323 << "ASCII\n"
3324 << "DATASET UNSTRUCTURED_GRID\n"
3325 << "POINTS " << n_vertices << " double\n";
3326
3327 // actually write the vertices.
3328 for (const auto &v : vertices)
3329 {
3330 out << v;
3331 for (unsigned int d = spacedim + 1; d <= 3; ++d)
3332 out << " 0"; // fill with zeroes
3333 out << '\n';
3334 }
3335
3336 const auto faces = vtk_flags.output_only_relevant ?
3337 get_relevant_face_iterators(tria) :
3338 get_boundary_face_iterators(tria);
3339 const auto edges = vtk_flags.output_only_relevant ?
3340 get_relevant_edge_iterators(tria) :
3341 get_boundary_edge_iterators(tria);
3342
3344 vtk_flags.output_cells || (dim >= 2 && vtk_flags.output_faces) ||
3345 (dim >= 3 && vtk_flags.output_edges),
3346 ExcMessage(
3347 "At least one of the flags (output_cells, output_faces, output_edges) has to be enabled!"));
3348
3349 // Write cells preamble
3350 const int n_cells = (vtk_flags.output_cells ? tria.n_active_cells() : 0) +
3351 (vtk_flags.output_faces ? faces.size() : 0) +
3352 (vtk_flags.output_edges ? edges.size() : 0);
3353
3354 // VTK now expects a number telling the total storage requirement to read all
3355 // cell connectivity information. The connectivity information is read cell by
3356 // cell, first specifying how many vertices are required to describe the cell,
3357 // and then specifying the index of every vertex. This means that for every
3358 // deal.II object type, we always need n_vertices + 1 integer per cell.
3359 // Compute the total number here.
3360 int cells_size = 0;
3361
3363 for (const auto &cell : tria.active_cell_iterators())
3364 cells_size += cell->n_vertices() + 1;
3365
3367 for (const auto &face : faces)
3368 cells_size += face->n_vertices() + 1;
3369
3371 for (const auto &edge : edges)
3372 cells_size += edge->n_vertices() + 1;
3373
3374 AssertThrow(cells_size > 0, ExcMessage("No cells given to be output!"));
3375
3376 out << "\nCELLS " << n_cells << ' ' << cells_size << '\n';
3377 /*
3378 * VTK cells:
3379 *
3380 * 1 VTK_VERTEX
3381 * 3 VTK_LINE
3382 * 5 VTK_TRIANGLE
3383 * 9 VTK_QUAD
3384 * 10 VTK_TETRA
3385 * 14 VTK_PYRAMID
3386 * 13 VTK_WEDGE
3387 * 12 VTK_HEXAHEDRON
3388 *
3389 * see also: https://vtk.org/wp-content/uploads/2015/04/file-formats.pdf
3390 */
3391 static const std::array<int, 8> deal_to_vtk_cell_type = {
3392 {1, 3, 5, 9, 10, 14, 13, 12}};
3393 static const std::array<unsigned int, 8> vtk_to_deal_hypercube = {
3394 {0, 1, 3, 2, 4, 5, 7, 6}};
3395
3396 // write cells.
3398 for (const auto &cell : tria.active_cell_iterators())
3399 {
3400 out << cell->n_vertices();
3401 for (const unsigned int i : cell->vertex_indices())
3402 {
3403 out << ' ';
3404 const auto reference_cell = cell->reference_cell();
3405
3406 if ((reference_cell == ReferenceCells::Vertex) ||
3407 (reference_cell == ReferenceCells::Line) ||
3408 (reference_cell == ReferenceCells::Quadrilateral) ||
3409 (reference_cell == ReferenceCells::Hexahedron))
3410 out << cell->vertex_index(vtk_to_deal_hypercube[i]);
3411 else if ((reference_cell == ReferenceCells::Triangle) ||
3412 (reference_cell == ReferenceCells::Tetrahedron) ||
3413 (reference_cell == ReferenceCells::Wedge))
3414 out << cell->vertex_index(i);
3415 else if (reference_cell == ReferenceCells::Pyramid)
3416 {
3417 static const std::array<unsigned int, 5> permutation_table{
3418 {0, 1, 3, 2, 4}};
3419 out << cell->vertex_index(permutation_table[i]);
3420 }
3421 else
3423 }
3424 out << '\n';
3425 }
3427 for (const auto &face : faces)
3428 {
3429 out << face->n_vertices();
3430 for (const unsigned int i : face->vertex_indices())
3431 {
3432 out << ' '
3433 << face->vertex_index(GeometryInfo<dim>::vertices_per_face ==
3434 face->n_vertices() ?
3435 vtk_to_deal_hypercube[i] :
3436 i);
3437 }
3438 out << '\n';
3439 }
3441 for (const auto &edge : edges)
3442 {
3443 out << 2;
3444 for (const unsigned int i : edge->vertex_indices())
3445 out << ' ' << edge->vertex_index(i);
3446 out << '\n';
3447 }
3448
3449 // write cell types
3450 out << "\nCELL_TYPES " << n_cells << '\n';
3452 {
3453 for (const auto &cell : tria.active_cell_iterators())
3454 out << deal_to_vtk_cell_type[static_cast<int>(cell->reference_cell())]
3455 << ' ';
3456 out << '\n';
3457 }
3459 {
3460 for (const auto &face : faces)
3461 out << deal_to_vtk_cell_type[static_cast<int>(face->reference_cell())]
3462 << ' ';
3463 out << '\n';
3464 }
3466 {
3467 for (const auto &edge : edges)
3468 out << deal_to_vtk_cell_type[static_cast<int>(edge->reference_cell())]
3469 << ' ';
3470 }
3471 out << "\n\nCELL_DATA " << n_cells << '\n'
3472 << "SCALARS MaterialID int 1\n"
3473 << "LOOKUP_TABLE default\n";
3474
3475 // Now material id and boundary id
3477 {
3478 for (const auto &cell : tria.active_cell_iterators())
3479 {
3480 out << static_cast<std::make_signed_t<types::material_id>>(
3481 cell->material_id())
3482 << ' ';
3483 }
3484 out << '\n';
3485 }
3487 {
3488 for (const auto &face : faces)
3489 {
3490 out << static_cast<std::make_signed_t<types::boundary_id>>(
3491 face->boundary_id())
3492 << ' ';
3493 }
3494 out << '\n';
3495 }
3497 {
3498 for (const auto &edge : edges)
3499 {
3500 out << static_cast<std::make_signed_t<types::boundary_id>>(
3501 edge->boundary_id())
3502 << ' ';
3503 }
3504 }
3505
3506 out << "\n\nSCALARS ManifoldID int 1\n"
3507 << "LOOKUP_TABLE default\n";
3508
3509 // Now manifold id
3511 {
3512 for (const auto &cell : tria.active_cell_iterators())
3513 {
3514 out << static_cast<std::make_signed_t<types::manifold_id>>(
3515 cell->manifold_id())
3516 << ' ';
3517 }
3518 out << '\n';
3519 }
3521 {
3522 for (const auto &face : faces)
3523 {
3524 out << static_cast<std::make_signed_t<types::manifold_id>>(
3525 face->manifold_id())
3526 << ' ';
3527 }
3528 out << '\n';
3529 }
3531 {
3532 for (const auto &edge : edges)
3533 {
3534 out << static_cast<std::make_signed_t<types::manifold_id>>(
3535 edge->manifold_id())
3536 << ' ';
3537 }
3538 out << '\n';
3539 }
3540
3541 out.flush();
3542
3543 AssertThrow(out.fail() == false, ExcIO());
3544}
3545
3546
3547
3548template <int dim, int spacedim>
3549void
3551 std::ostream &out) const
3552{
3553 AssertThrow(out.fail() == false, ExcIO());
3554
3555 // convert the cells of the triangulation into a set of patches
3556 // and then have them output. since there is no data attached to
3557 // the geometry, we also do not have to provide any names, identifying
3558 // information, etc.
3559 std::vector<DataOutBase::Patch<dim, spacedim>> patches;
3560 patches.reserve(tria.n_active_cells());
3561 generate_triangulation_patches(patches, tria.begin_active(), tria.end());
3562
3565 patches,
3566 triangulation_patch_data_names(),
3567 std::vector<
3568 std::tuple<unsigned int,
3569 unsigned int,
3570 std::string,
3572 vtu_flags,
3573 out);
3575 {
3576 out << " </UnstructuredGrid>\n";
3577 out << "<dealiiData encoding=\"base64\">";
3578 std::stringstream outstring;
3579 boost::archive::binary_oarchive ia(outstring);
3580 tria.save(ia, 0);
3581 const auto compressed = Utilities::compress(outstring.str());
3582 out << Utilities::encode_base64({compressed.begin(), compressed.end()});
3583 out << "\n</dealiiData>\n";
3584 out << "</VTKFile>\n";
3585 }
3586 else
3588
3589 out << std::flush;
3590 AssertThrow(out.fail() == false, ExcIO());
3591}
3592
3593
3594
3595template <int dim, int spacedim>
3596void
3598 const Triangulation<dim, spacedim> &tria,
3599 const std::string &filename_without_extension,
3600 const bool view_levels,
3601 const bool include_artificial) const
3602{
3603 std::vector<DataOutBase::Patch<dim, spacedim>> patches;
3604 const unsigned int n_datasets = 4;
3605 std::vector<std::string> data_names;
3606 data_names.emplace_back("level");
3607 data_names.emplace_back("subdomain");
3608 data_names.emplace_back("level_subdomain");
3609 data_names.emplace_back("proc_writing");
3610
3611 const auto &reference_cells = tria.get_reference_cells();
3612
3613 AssertDimension(reference_cells.size(), 1);
3614
3615 const auto &reference_cell = reference_cells[0];
3616
3617 const unsigned int n_q_points = reference_cell.n_vertices();
3618
3619 for (const auto &cell : tria.cell_iterators())
3620 {
3621 if (!view_levels)
3622 {
3623 if (cell->has_children())
3624 continue;
3625 if (!include_artificial &&
3626 cell->subdomain_id() == numbers::artificial_subdomain_id)
3627 continue;
3628 }
3629 else if (!include_artificial)
3630 {
3631 if (cell->has_children() &&
3632 cell->level_subdomain_id() == numbers::artificial_subdomain_id)
3633 continue;
3634 else if (cell->is_active() &&
3635 cell->level_subdomain_id() ==
3637 cell->subdomain_id() == numbers::artificial_subdomain_id)
3638 continue;
3639 }
3640
3642 patch.data.reinit(n_datasets, n_q_points);
3643 patch.points_are_available = false;
3644 patch.reference_cell = reference_cell;
3645
3646 for (unsigned int vertex = 0; vertex < n_q_points; ++vertex)
3647 {
3648 patch.vertices[vertex] = cell->vertex(vertex);
3649 patch.data(0, vertex) = cell->level();
3650 if (cell->is_active())
3651 patch.data(1, vertex) = static_cast<double>(
3652 static_cast<std::make_signed_t<types::subdomain_id>>(
3653 cell->subdomain_id()));
3654 else
3655 patch.data(1, vertex) = -1.0;
3656 patch.data(2, vertex) = static_cast<double>(
3657 static_cast<std::make_signed_t<types::subdomain_id>>(
3658 cell->level_subdomain_id()));
3659 patch.data(3, vertex) = tria.locally_owned_subdomain();
3660 }
3661
3662 for (auto f : reference_cell.face_indices())
3664 patches.push_back(patch);
3665 }
3666
3667 // only create .pvtu file if running in parallel
3668 // if not, just create a .vtu file with no reference
3669 // to the processor number
3670 std::string new_file = filename_without_extension + ".vtu";
3672 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(&tria))
3673 {
3674 new_file = filename_without_extension + ".proc" +
3675 Utilities::int_to_string(tr->locally_owned_subdomain(), 4) +
3676 ".vtu";
3677
3678 // create .pvtu record
3679 if (tr->locally_owned_subdomain() == 0)
3680 {
3681 std::vector<std::string> filenames;
3682
3683 // .pvtu needs to reference the files without a relative path because
3684 // it will be written in the same directory. For this, remove any
3685 // paths from filename.
3686 std::size_t pos = filename_without_extension.find_last_of('/');
3687 if (pos == std::string::npos)
3688 pos = 0;
3689 else
3690 pos += 1;
3691 const unsigned int n_procs =
3692 Utilities::MPI::n_mpi_processes(tr->get_mpi_communicator());
3693 filenames.reserve(n_procs);
3694 for (unsigned int i = 0; i < n_procs; ++i)
3695 filenames.push_back(filename_without_extension.substr(pos) +
3696 ".proc" + Utilities::int_to_string(i, 4) +
3697 ".vtu");
3698
3699 const std::string pvtu_filename =
3700 (filename_without_extension + ".pvtu");
3701 std::ofstream pvtu_output(pvtu_filename);
3702
3703 DataOut<dim, spacedim> data_out;
3704 data_out.attach_triangulation(*tr);
3705
3706 // We need a dummy vector with the names of the data values in the
3707 // .vtu files in order that the .pvtu contains reference these values
3708 const Vector<float> dummy_vector(tr->n_active_cells());
3709 data_out.add_data_vector(dummy_vector, "level");
3710 data_out.add_data_vector(dummy_vector, "subdomain");
3711 data_out.add_data_vector(dummy_vector, "level_subdomain");
3712 data_out.add_data_vector(dummy_vector, "proc_writing");
3713
3714 data_out.build_patches();
3715
3716 data_out.write_pvtu_record(pvtu_output, filenames);
3717 }
3718 }
3719
3720 std::ofstream out(new_file);
3721 std::vector<
3722 std::tuple<unsigned int,
3723 unsigned int,
3724 std::string,
3726 vector_data_ranges;
3728 patches, data_names, vector_data_ranges, vtu_flags, out);
3729}
3730
3731
3732
3733unsigned int
3735{
3736 return 0;
3737}
3738
3739unsigned int
3741{
3742 return 0;
3743}
3744
3745
3746unsigned int
3748{
3749 return 0;
3750}
3751
3752unsigned int
3754{
3755 return 0;
3756}
3757
3758unsigned int
3760{
3761 return 0;
3762}
3763
3764unsigned int
3766{
3767 return 0;
3768}
3769
3770unsigned int
3772{
3773 return 0;
3774}
3775
3776unsigned int
3778{
3779 return 0;
3780}
3781
3782
3783
3784template <int dim, int spacedim>
3785unsigned int
3787{
3789 unsigned int n_faces = 0;
3790
3791 for (const auto &face : tria.active_face_iterators())
3792 if ((face->at_boundary()) && (face->boundary_id() != 0))
3793 ++n_faces;
3794
3795 return n_faces;
3796}
3797
3798
3799
3800template <int dim, int spacedim>
3801unsigned int
3803{
3804 // save the user flags for lines so
3805 // we can use these flags to track
3806 // which ones we've already counted
3807 std::vector<bool> line_flags;
3808 const_cast<::Triangulation<dim, spacedim> &>(tria).save_user_flags_line(
3809 line_flags);
3810 const_cast<::Triangulation<dim, spacedim> &>(tria)
3811 .clear_user_flags_line();
3812
3813 unsigned int n_lines = 0;
3814
3815 for (const auto &cell : tria.active_cell_iterators())
3816 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
3817 if (cell->line(l)->at_boundary() && (cell->line(l)->boundary_id() != 0) &&
3818 (cell->line(l)->user_flag_set() == false))
3819 {
3820 ++n_lines;
3821 cell->line(l)->set_user_flag();
3822 }
3823
3824 // at the end, restore the user
3825 // flags for the lines
3826 const_cast<::Triangulation<dim, spacedim> &>(tria).load_user_flags_line(
3827 line_flags);
3828
3829 return n_lines;
3830}
3831
3832
3833
3834unsigned int
3836 const unsigned int next_element_index,
3837 std::ostream &) const
3838{
3839 return next_element_index;
3840}
3841
3842
3843unsigned int
3845 const unsigned int next_element_index,
3846 std::ostream &) const
3847{
3848 return next_element_index;
3849}
3850
3851unsigned int
3853 const unsigned int next_element_index,
3854 std::ostream &) const
3855{
3856 return next_element_index;
3857}
3858
3859
3860unsigned int
3862 const unsigned int next_element_index,
3863 std::ostream &) const
3864{
3865 return next_element_index;
3866}
3867
3868unsigned int
3870 const unsigned int next_element_index,
3871 std::ostream &) const
3872{
3873 return next_element_index;
3874}
3875
3876
3877unsigned int
3879 const unsigned int next_element_index,
3880 std::ostream &) const
3881{
3882 return next_element_index;
3883}
3884
3885
3886unsigned int
3888 const unsigned int next_element_index,
3889 std::ostream &) const
3890{
3891 return next_element_index;
3892}
3893
3894unsigned int
3896 const unsigned int next_element_index,
3897 std::ostream &) const
3898{
3899 return next_element_index;
3900}
3901
3902
3903
3904template <int dim, int spacedim>
3905unsigned int
3907 const unsigned int next_element_index,
3908 std::ostream &out) const
3909{
3910 unsigned int current_element_index = next_element_index;
3911
3912 for (const auto &face : tria.active_face_iterators())
3913 if (face->at_boundary() && (face->boundary_id() != 0))
3914 {
3915 out << current_element_index << ' '
3916 << face->reference_cell().gmsh_element_type() << ' ';
3917 out << static_cast<unsigned int>(face->boundary_id()) << ' '
3918 << static_cast<unsigned int>(face->boundary_id()) << ' '
3919 << face->n_vertices();
3920 // note: vertex numbers are 1-base
3921 for (const unsigned int vertex : face->vertex_indices())
3922 {
3923 if (face->reference_cell() == ReferenceCells::Quadrilateral)
3924 out << ' '
3925 << face->vertex_index(
3927 1;
3928 else if ((face->reference_cell() == ReferenceCells::Triangle) ||
3929 (face->reference_cell() == ReferenceCells::Line))
3930 out << ' ' << face->vertex_index(vertex) + 1;
3931 else
3933 }
3934 out << '\n';
3935
3936 ++current_element_index;
3937 }
3938 return current_element_index;
3939}
3940
3941
3942
3943template <int dim, int spacedim>
3944unsigned int
3946 const unsigned int next_element_index,
3947 std::ostream &out) const
3948{
3949 unsigned int current_element_index = next_element_index;
3950 // save the user flags for lines so
3951 // we can use these flags to track
3952 // which ones we've already taken
3953 // care of
3954 std::vector<bool> line_flags;
3955 const_cast<::Triangulation<dim, spacedim> &>(tria).save_user_flags_line(
3956 line_flags);
3957 const_cast<::Triangulation<dim, spacedim> &>(tria)
3958 .clear_user_flags_line();
3959
3960 for (const auto &cell : tria.active_cell_iterators())
3961 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
3962 if (cell->line(l)->at_boundary() && (cell->line(l)->boundary_id() != 0) &&
3963 (cell->line(l)->user_flag_set() == false))
3964 {
3965 out << next_element_index << ' '
3967 out << static_cast<unsigned int>(cell->line(l)->boundary_id()) << ' '
3968 << static_cast<unsigned int>(cell->line(l)->boundary_id())
3969 << " 2 "; // two vertex indices to follow
3970 // note: vertex numbers are 1-base
3971 for (unsigned int vertex = 0; vertex < 2; ++vertex)
3972 out << ' '
3973 << cell->line(l)->vertex_index(
3975 1;
3976 out << '\n';
3977
3978 // move on to the next line
3979 // but mark the current one
3980 // as taken care of
3981 ++current_element_index;
3982 cell->line(l)->set_user_flag();
3983 }
3984
3985 // at the end, restore the user
3986 // flags for the lines
3987 const_cast<::Triangulation<dim, spacedim> &>(tria).load_user_flags_line(
3988 line_flags);
3989
3990 return current_element_index;
3991}
3992
3993
3994
3995unsigned int
3997 const unsigned int next_element_index,
3998 std::ostream &) const
3999{
4000 return next_element_index;
4001}
4002
4003unsigned int
4005 const unsigned int next_element_index,
4006 std::ostream &) const
4007{
4008 return next_element_index;
4009}
4010
4011unsigned int
4013 const unsigned int next_element_index,
4014 std::ostream &) const
4015{
4016 return next_element_index;
4017}
4018
4019unsigned int
4021 const unsigned int next_element_index,
4022 std::ostream &) const
4023{
4024 return next_element_index;
4025}
4026
4027unsigned int
4029 const unsigned int next_element_index,
4030 std::ostream &) const
4031{
4032 return next_element_index;
4033}
4034
4035
4036unsigned int
4038 const unsigned int next_element_index,
4039 std::ostream &) const
4040{
4041 return next_element_index;
4042}
4043
4044
4045unsigned int
4047 const unsigned int next_element_index,
4048 std::ostream &) const
4049{
4050 return next_element_index;
4051}
4052
4053unsigned int
4055 const unsigned int next_element_index,
4056 std::ostream &) const
4057{
4058 return next_element_index;
4059}
4060
4061
4062
4063template <int dim, int spacedim>
4064unsigned int
4066 const unsigned int next_element_index,
4067 std::ostream &out) const
4068{
4069 unsigned int current_element_index = next_element_index;
4071
4072 for (const auto &face : tria.active_face_iterators())
4073 if (face->at_boundary() && (face->boundary_id() != 0))
4074 {
4075 out << current_element_index << " "
4076 << static_cast<unsigned int>(face->boundary_id()) << " ";
4077 switch (dim)
4078 {
4079 case 2:
4080 out << "line ";
4081 break;
4082 case 3:
4083 out << "quad ";
4084 break;
4085 default:
4087 }
4088 // note: vertex numbers are 1-base
4089 for (unsigned int vertex = 0;
4090 vertex < GeometryInfo<dim>::vertices_per_face;
4091 ++vertex)
4092 out << face->vertex_index(
4094 1
4095 << ' ';
4096 out << '\n';
4097
4098 ++current_element_index;
4099 }
4100 return current_element_index;
4101}
4102
4103
4104
4105template <int dim, int spacedim>
4106unsigned int
4108 const unsigned int next_element_index,
4109 std::ostream &out) const
4110{
4111 unsigned int current_element_index = next_element_index;
4112 // save the user flags for lines so
4113 // we can use these flags to track
4114 // which ones we've already taken
4115 // care of
4116 std::vector<bool> line_flags;
4117 const_cast<::Triangulation<dim, spacedim> &>(tria).save_user_flags_line(
4118 line_flags);
4119 const_cast<::Triangulation<dim, spacedim> &>(tria)
4120 .clear_user_flags_line();
4121
4122 for (const auto &cell : tria.active_cell_iterators())
4123 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
4124 if (cell->line(l)->at_boundary() && (cell->line(l)->boundary_id() != 0) &&
4125 (cell->line(l)->user_flag_set() == false))
4126 {
4127 out << current_element_index << " "
4128 << static_cast<unsigned int>(cell->line(l)->boundary_id())
4129 << " line ";
4130 // note: vertex numbers in ucd format are 1-base
4131 for (unsigned int vertex = 0; vertex < 2; ++vertex)
4132 out << cell->line(l)->vertex_index(
4134 1
4135 << ' ';
4136 out << '\n';
4137
4138 // move on to the next line
4139 // but mark the current one
4140 // as taken care of
4141 ++current_element_index;
4142 cell->line(l)->set_user_flag();
4143 }
4144
4145 // at the end, restore the user
4146 // flags for the lines
4147 const_cast<::Triangulation<dim, spacedim> &>(tria).load_user_flags_line(
4148 line_flags);
4149 return current_element_index;
4150}
4151
4152
4153
4154namespace internal
4155{
4156 namespace
4157 {
4166 template <int spacedim>
4167 void
4168 remove_collinear_points(std::vector<Point<spacedim>> &points)
4169 {
4170 while (points.size() > 2)
4171 {
4172 Tensor<1, spacedim> first_difference = points[1] - points[0];
4173 first_difference /= first_difference.norm();
4174 Tensor<1, spacedim> second_difference = points[2] - points[1];
4175 second_difference /= second_difference.norm();
4176 // If the three points are collinear then remove the middle one.
4177 if ((first_difference - second_difference).norm() < 1e-10)
4178 points.erase(points.begin() + 1);
4179 else
4180 break;
4181 }
4182 }
4183
4184
4185
4186 template <int spacedim>
4187 void
4188 write_gnuplot(const ::Triangulation<1, spacedim> &tria,
4189 std::ostream &out,
4190 const Mapping<1, spacedim> *,
4191 const GridOutFlags::Gnuplot &gnuplot_flags)
4192 {
4193 AssertThrow(out.fail() == false, ExcIO());
4194
4195 for (const auto &cell : tria.active_cell_iterators())
4196 {
4197 if (gnuplot_flags.write_cell_numbers)
4198 out << "# cell " << cell << '\n';
4199
4200 out << cell->vertex(0) << ' ' << cell->level() << ' '
4201 << cell->material_id() << '\n'
4202 << cell->vertex(1) << ' ' << cell->level() << ' '
4203 << cell->material_id() << '\n'
4204 << "\n\n";
4205 }
4206
4207 // make sure everything now gets to
4208 // disk
4209 out.flush();
4210
4211 AssertThrow(out.fail() == false, ExcIO());
4212 }
4213
4214
4215
4216 template <int spacedim>
4217 void
4218 write_gnuplot(const ::Triangulation<2, spacedim> &tria,
4219 std::ostream &out,
4220 const Mapping<2, spacedim> *mapping,
4221 const GridOutFlags::Gnuplot &gnuplot_flags)
4222 {
4223 AssertThrow(out.fail() == false, ExcIO());
4224
4225 const int dim = 2;
4226
4227 const unsigned int n_additional_points =
4228 gnuplot_flags.n_extra_curved_line_points;
4229 const unsigned int n_points = 2 + n_additional_points;
4230
4231 // If we need to plot curved lines then generate a quadrature formula to
4232 // place points via the mapping
4233 Quadrature<dim> q_projector;
4234 std::vector<Point<dim - 1>> boundary_points;
4235 if (mapping != nullptr)
4236 {
4237 boundary_points.resize(n_points);
4238 boundary_points[0][0] = 0;
4239 boundary_points[n_points - 1][0] = 1;
4240 for (unsigned int i = 1; i < n_points - 1; ++i)
4241 boundary_points[i][0] = 1. * i / (n_points - 1);
4242
4243 const std::vector<double> dummy_weights(n_points, 1. / n_points);
4244 const Quadrature<dim - 1> quadrature(boundary_points, dummy_weights);
4245
4247 ReferenceCells::get_hypercube<dim>(), quadrature);
4248 }
4249
4250 static constexpr std::array<unsigned int, 8> local_vertex_numbering = {
4251 {0, 1, 5, 4, 2, 3, 7, 6}};
4252 for (const auto &cell : tria.active_cell_iterators())
4253 {
4254 if (gnuplot_flags.write_cell_numbers)
4255 out << "# cell " << cell << '\n';
4256
4257 if (mapping == nullptr ||
4258 (dim == spacedim ?
4259 (!cell->at_boundary() && !gnuplot_flags.curved_inner_cells) :
4260 // ignore checking for boundary or interior cells in the codim
4261 // 1 case: 'or false' is a no-op
4262 false))
4263 {
4264 // write out the four sides of this cell by putting the four
4265 // points (+ the initial point again) in a row and lifting the
4266 // drawing pencil at the end
4267 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
4268 out << cell->vertex(dim == 3 ?
4269 local_vertex_numbering[i] :
4270 GeometryInfo<dim>::ucd_to_deal[i])
4271 << ' ' << cell->level() << ' ' << cell->material_id()
4272 << '\n';
4273 out << cell->vertex(0) << ' ' << cell->level() << ' '
4274 << cell->material_id() << '\n'
4275 << '\n' // double new line for gnuplot 3d plots
4276 << '\n';
4277 }
4278 else
4279 // cell is at boundary and we are to treat curved boundaries. so
4280 // loop over all faces and draw them as small pieces of lines
4281 {
4282 for (const unsigned int face_no :
4283 GeometryInfo<dim>::face_indices())
4284 {
4285 const typename ::Triangulation<dim,
4286 spacedim>::face_iterator
4287 face = cell->face(face_no);
4288 if (dim != spacedim || face->at_boundary() ||
4289 gnuplot_flags.curved_inner_cells)
4290 {
4291 // Save the points on each face to a vector and then try
4292 // to remove collinear points that won't show up in the
4293 // generated plot.
4294 std::vector<Point<spacedim>> line_points;
4295 // compute offset of quadrature points within set of
4296 // projected points
4297 const auto offset =
4299 cell->reference_cell(),
4300 face_no,
4301 cell->combined_face_orientation(face_no),
4302 n_points);
4303 line_points.reserve(n_points);
4304 for (unsigned int i = 0; i < n_points; ++i)
4305 line_points.push_back(
4307 cell, q_projector.point(offset + i)));
4308 internal::remove_collinear_points(line_points);
4309
4310 for (const Point<spacedim> &point : line_points)
4311 out << point << ' ' << cell->level() << ' '
4312 << cell->material_id() << '\n';
4313
4314 out << '\n' << '\n';
4315 }
4316 else
4317 {
4318 // if, however, the face is not at the boundary and we
4319 // don't want to curve anything, then draw it as usual
4320 out << face->vertex(0) << ' ' << cell->level() << ' '
4321 << cell->material_id() << '\n'
4322 << face->vertex(1) << ' ' << cell->level() << ' '
4323 << cell->material_id() << '\n'
4324 << '\n'
4325 << '\n';
4326 }
4327 }
4328 }
4329 }
4330
4331 // make sure everything now gets to disk
4332 out.flush();
4333
4334 AssertThrow(out.fail() == false, ExcIO());
4335 }
4336
4337
4338
4339 template <int spacedim>
4340 void
4341 write_gnuplot(const ::Triangulation<3, spacedim> &tria,
4342 std::ostream &out,
4343 const Mapping<3, spacedim> *mapping,
4344 const GridOutFlags::Gnuplot &gnuplot_flags)
4345 {
4346 AssertThrow(out.fail() == false, ExcIO());
4347
4348 const int dim = 3;
4349
4350 const unsigned int n_additional_points =
4351 gnuplot_flags.n_extra_curved_line_points;
4352 const unsigned int n_points = 2 + n_additional_points;
4353
4354 // If we need to plot curved lines then generate a quadrature formula to
4355 // place points via the mapping
4356 std::unique_ptr<Quadrature<dim>> q_projector;
4357 std::vector<Point<1>> boundary_points;
4358 if (mapping != nullptr)
4359 {
4360 boundary_points.resize(n_points);
4361 boundary_points[0][0] = 0;
4362 boundary_points[n_points - 1][0] = 1;
4363 for (unsigned int i = 1; i < n_points - 1; ++i)
4364 boundary_points[i][0] = 1. * i / (n_points - 1);
4365
4366 const std::vector<double> dummy_weights(n_points, 1. / n_points);
4367 const Quadrature<1> quadrature1d(boundary_points, dummy_weights);
4368
4369 // tensor product of points, only one copy
4370 const QIterated<dim - 1> quadrature(quadrature1d, 1);
4371 q_projector = std::make_unique<Quadrature<dim>>(
4373 ReferenceCells::get_hypercube<dim>(), quadrature));
4374 }
4375
4376 for (const auto &cell : tria.active_cell_iterators())
4377 {
4378 if (gnuplot_flags.write_cell_numbers)
4379 out << "# cell " << cell << '\n';
4380
4381 if (mapping == nullptr || n_points == 2 ||
4382 (!cell->has_boundary_lines() &&
4383 !gnuplot_flags.curved_inner_cells))
4384 {
4385 if (cell->reference_cell() == ReferenceCells::Hexahedron)
4386 {
4387 // front face
4388 out << cell->vertex(0) << ' ' << cell->level() << ' '
4389 << cell->material_id() << '\n'
4390 << cell->vertex(1) << ' ' << cell->level() << ' '
4391 << cell->material_id() << '\n'
4392 << cell->vertex(5) << ' ' << cell->level() << ' '
4393 << cell->material_id() << '\n'
4394 << cell->vertex(4) << ' ' << cell->level() << ' '
4395 << cell->material_id() << '\n'
4396 << cell->vertex(0) << ' ' << cell->level() << ' '
4397 << cell->material_id() << '\n'
4398 << '\n';
4399 // back face
4400 out << cell->vertex(2) << ' ' << cell->level() << ' '
4401 << cell->material_id() << '\n'
4402 << cell->vertex(3) << ' ' << cell->level() << ' '
4403 << cell->material_id() << '\n'
4404 << cell->vertex(7) << ' ' << cell->level() << ' '
4405 << cell->material_id() << '\n'
4406 << cell->vertex(6) << ' ' << cell->level() << ' '
4407 << cell->material_id() << '\n'
4408 << cell->vertex(2) << ' ' << cell->level() << ' '
4409 << cell->material_id() << '\n'
4410 << '\n';
4411
4412 // now for the four connecting lines
4413 out << cell->vertex(0) << ' ' << cell->level() << ' '
4414 << cell->material_id() << '\n'
4415 << cell->vertex(2) << ' ' << cell->level() << ' '
4416 << cell->material_id() << '\n'
4417 << '\n';
4418 out << cell->vertex(1) << ' ' << cell->level() << ' '
4419 << cell->material_id() << '\n'
4420 << cell->vertex(3) << ' ' << cell->level() << ' '
4421 << cell->material_id() << '\n'
4422 << '\n';
4423 out << cell->vertex(5) << ' ' << cell->level() << ' '
4424 << cell->material_id() << '\n'
4425 << cell->vertex(7) << ' ' << cell->level() << ' '
4426 << cell->material_id() << '\n'
4427 << '\n';
4428 out << cell->vertex(4) << ' ' << cell->level() << ' '
4429 << cell->material_id() << '\n'
4430 << cell->vertex(6) << ' ' << cell->level() << ' '
4431 << cell->material_id() << '\n'
4432 << '\n';
4433 }
4434 else if (cell->reference_cell() == ReferenceCells::Tetrahedron)
4435 {
4436 // Draw the tetrahedron as a two collections of lines.
4437 for (const unsigned int v : {0, 1, 2, 0, 3, 2})
4438 {
4439 out << cell->vertex(v) << ' ' << cell->level() << ' '
4440 << cell->material_id() << '\n';
4441 }
4442 out << '\n'; // end of first line
4443
4444 for (const unsigned int v : {3, 1})
4445 {
4446 out << cell->vertex(v) << ' ' << cell->level() << ' '
4447 << cell->material_id() << '\n';
4448 }
4449 out << '\n'; // end of second line
4450 }
4451 else if (cell->reference_cell() == ReferenceCells::Wedge)
4452 {
4453 // Draw the wedge as a collection of three
4454 // lines. The first one wraps around the base,
4455 // goes up to the top, and wraps around that. The
4456 // second and third are just individual lines
4457 // going from base to top.
4458 for (const unsigned int v : {0, 1, 2, 0, 3, 4, 5, 3})
4459 {
4460 out << cell->vertex(v) << ' ' << cell->level() << ' '
4461 << cell->material_id() << '\n';
4462 }
4463 out << '\n'; // end of first line
4464
4465 for (const unsigned int v : {1, 4})
4466 {
4467 out << cell->vertex(v) << ' ' << cell->level() << ' '
4468 << cell->material_id() << '\n';
4469 }
4470 out << '\n'; // end of second line
4471
4472 for (const unsigned int v : {2, 5})
4473 {
4474 out << cell->vertex(v) << ' ' << cell->level() << ' '
4475 << cell->material_id() << '\n';
4476 }
4477 out << '\n'; // end of third line
4478 }
4479 else if (cell->reference_cell() == ReferenceCells::Pyramid)
4480 {
4481 // Draw the pyramid as a collections of two lines.
4482 for (const unsigned int v : {0, 1, 3, 2, 0, 4, 1})
4483 {
4484 out << cell->vertex(v) << ' ' << cell->level() << ' '
4485 << cell->material_id() << '\n';
4486 }
4487 out << '\n'; // end of first line
4488
4489 for (const unsigned int v : {2, 4, 3})
4490 {
4491 out << cell->vertex(v) << ' ' << cell->level() << ' '
4492 << cell->material_id() << '\n';
4493 }
4494 out << '\n'; // end of second line
4495 }
4496 else
4498 }
4499 else // need to handle curved boundaries
4500 {
4501 Assert(cell->reference_cell() == ReferenceCells::Hexahedron,
4503 for (const unsigned int face_no :
4504 GeometryInfo<dim>::face_indices())
4505 {
4506 const typename ::Triangulation<dim,
4507 spacedim>::face_iterator
4508 face = cell->face(face_no);
4509
4510 if (face->at_boundary() &&
4511 gnuplot_flags.write_additional_boundary_lines)
4512 {
4513 const auto offset =
4515 cell->reference_cell(),
4516 face_no,
4517 cell->combined_face_orientation(face_no),
4518 n_points * n_points);
4519 for (unsigned int i = 0; i < n_points - 1; ++i)
4520 for (unsigned int j = 0; j < n_points - 1; ++j)
4521 {
4522 const Point<spacedim> p0 =
4524 cell,
4525 q_projector->point(offset + i * n_points + j));
4526 out << p0 << ' ' << cell->level() << ' '
4527 << cell->material_id() << '\n';
4528 out << (mapping->transform_unit_to_real_cell(
4529 cell,
4530 q_projector->point(
4531 offset + (i + 1) * n_points + j)))
4532 << ' ' << cell->level() << ' '
4533 << cell->material_id() << '\n';
4534 out << (mapping->transform_unit_to_real_cell(
4535 cell,
4536 q_projector->point(
4537 offset + (i + 1) * n_points + j + 1)))
4538 << ' ' << cell->level() << ' '
4539 << cell->material_id() << '\n';
4540 out << (mapping->transform_unit_to_real_cell(
4541 cell,
4542 q_projector->point(offset + i * n_points +
4543 j + 1)))
4544 << ' ' << cell->level() << ' '
4545 << cell->material_id() << '\n';
4546 // and the first point again
4547 out << p0 << ' ' << cell->level() << ' '
4548 << cell->material_id() << '\n';
4549 out << '\n' << '\n';
4550 }
4551 }
4552 else
4553 {
4554 for (unsigned int l = 0;
4555 l < GeometryInfo<dim>::lines_per_face;
4556 ++l)
4557 {
4558 const typename ::Triangulation<dim, spacedim>::
4559 line_iterator line = face->line(l);
4560
4561 const Point<spacedim> &v0 = line->vertex(0),
4562 &v1 = line->vertex(1);
4563 if (line->at_boundary() ||
4564 gnuplot_flags.curved_inner_cells)
4565 {
4566 // Save the points on each face to a vector and
4567 // then try to remove collinear points that won't
4568 // show up in the generated plot.
4569 std::vector<Point<spacedim>> line_points;
4570 // transform_real_to_unit_cell could be replaced
4571 // by using QProjector<dim>::project_to_line
4572 // which is not yet implemented
4573 const Point<spacedim>
4574 u0 = mapping->transform_real_to_unit_cell(cell,
4575 v0),
4576 u1 = mapping->transform_real_to_unit_cell(cell,
4577 v1);
4578 line_points.reserve(n_points);
4579 for (unsigned int i = 0; i < n_points; ++i)
4580 line_points.push_back(
4582 cell,
4583 (1 - boundary_points[i][0]) * u0 +
4584 boundary_points[i][0] * u1));
4585 internal::remove_collinear_points(line_points);
4586 for (const Point<spacedim> &point : line_points)
4587 out << point << ' ' << cell->level() << ' '
4588 << static_cast<unsigned int>(
4589 cell->material_id())
4590 << '\n';
4591 }
4592 else
4593 out << v0 << ' ' << cell->level() << ' '
4594 << cell->material_id() << '\n'
4595 << v1 << ' ' << cell->level() << ' '
4596 << cell->material_id() << '\n';
4597
4598 out << '\n' << '\n';
4599 }
4600 }
4601 }
4602 }
4603 }
4604
4605 // make sure everything now gets to disk
4606 out.flush();
4607
4608 AssertThrow(out.fail() == false, ExcIO());
4609 }
4610 } // namespace
4611} // namespace internal
4612
4613
4614
4615template <int dim, int spacedim>
4616void
4618 std::ostream &out,
4619 const Mapping<dim, spacedim> *mapping) const
4620{
4621 internal::write_gnuplot(tria, out, mapping, gnuplot_flags);
4622}
4623
4624
4625
4626namespace internal
4627{
4628 namespace
4629 {
4630 struct LineEntry
4631 {
4635 unsigned int level;
4636 LineEntry(const Point<2> &f,
4637 const Point<2> &s,
4638 const bool c,
4639 const unsigned int l)
4640 : first(f)
4641 , second(s)
4642 , colorize(c)
4643 , level(l)
4644 {}
4645 };
4646
4647
4648 void
4649 write_eps(const ::Triangulation<1> &,
4650 std::ostream &,
4651 const Mapping<1> *,
4652 const GridOutFlags::Eps<2> &,
4653 const GridOutFlags::Eps<3> &)
4654 {
4656 }
4657
4658 void
4659 write_eps(const ::Triangulation<1, 2> &,
4660 std::ostream &,
4661 const Mapping<1, 2> *,
4662 const GridOutFlags::Eps<2> &,
4663 const GridOutFlags::Eps<3> &)
4664 {
4666 }
4667
4668 void
4669 write_eps(const ::Triangulation<1, 3> &,
4670 std::ostream &,
4671 const Mapping<1, 3> *,
4672 const GridOutFlags::Eps<2> &,
4673 const GridOutFlags::Eps<3> &)
4674 {
4676 }
4677
4678 void
4679 write_eps(const ::Triangulation<2, 3> &,
4680 std::ostream &,
4681 const Mapping<2, 3> *,
4682 const GridOutFlags::Eps<2> &,
4683 const GridOutFlags::Eps<3> &)
4684 {
4686 }
4687
4688
4689
4690 template <int dim, int spacedim>
4691 void
4692 write_eps(const ::Triangulation<dim, spacedim> &tria,
4693 std::ostream &out,
4694 const Mapping<dim, spacedim> *mapping,
4695 const GridOutFlags::Eps<2> &eps_flags_2,
4696 const GridOutFlags::Eps<3> &eps_flags_3)
4697 {
4698 using LineList = std::list<LineEntry>;
4699
4700 // We should never get here in 1d since this function is overloaded for
4701 // all dim == 1 cases.
4702 Assert(dim == 2 || dim == 3, ExcInternalError());
4703
4704 // Copy, with an object slice, something containing the flags common to
4705 // all dimensions in order to avoid the recurring distinctions between
4706 // the different eps_flags present.
4707 const GridOutFlags::EpsFlagsBase eps_flags_base =
4708 dim == 2 ?
4709 static_cast<const GridOutFlags::EpsFlagsBase &>(eps_flags_2) :
4710 static_cast<const GridOutFlags::EpsFlagsBase &>(eps_flags_3);
4711
4712 AssertThrow(out.fail() == false, ExcIO());
4713 const unsigned int n_points = eps_flags_base.n_boundary_face_points;
4714
4715 // make up a list of lines by which
4716 // we will construct the triangulation
4717 //
4718 // this part unfortunately is a bit
4719 // dimension dependent, so we have to
4720 // treat every dimension different.
4721 // however, by directly producing
4722 // the lines to be printed, i.e. their
4723 // 2d images, we can later do the
4724 // actual output dimension independent
4725 // again
4726 LineList line_list;
4727
4728 switch (dim)
4729 {
4730 case 1:
4731 {
4733 break;
4734 }
4735
4736 case 2:
4737 {
4738 for (const auto &cell : tria.active_cell_iterators())
4739 for (const unsigned int line_no : cell->line_indices())
4740 {
4741 typename ::Triangulation<dim, spacedim>::line_iterator
4742 line = cell->line(line_no);
4743
4744 // first treat all
4745 // interior lines and
4746 // make up a list of
4747 // them. if curved
4748 // lines shall not be
4749 // supported (i.e. no
4750 // mapping is
4751 // provided), then also
4752 // treat all other
4753 // lines
4754 if (!line->has_children() &&
4755 (mapping == nullptr || !line->at_boundary()))
4756 // one would expect
4757 // make_pair(line->vertex(0),
4758 // line->vertex(1))
4759 // here, but that is
4760 // not dimension
4761 // independent, since
4762 // vertex(i) is
4763 // Point<dim>, but we
4764 // want a Point<2>.
4765 // in fact, whenever
4766 // we're here, the
4767 // vertex is a
4768 // Point<dim>, but
4769 // the compiler does
4770 // not know
4771 // this. hopefully,
4772 // the compiler will
4773 // optimize away this
4774 // little kludge
4775 line_list.emplace_back(
4776 Point<2>(line->vertex(0)[0], line->vertex(0)[1]),
4777 Point<2>(line->vertex(1)[0], line->vertex(1)[1]),
4778 line->user_flag_set(),
4779 cell->level());
4780 }
4781
4782 // next if we are to treat
4783 // curved boundaries
4784 // specially, then add lines
4785 // to the list consisting of
4786 // pieces of the boundary
4787 // lines
4788 if (mapping != nullptr)
4789 {
4790 // to do so, first
4791 // generate a sequence of
4792 // points on a face and
4793 // project them onto the
4794 // faces of a unit cell
4795 std::vector<Point<dim - 1>> boundary_points(n_points);
4796
4797 for (unsigned int i = 0; i < n_points; ++i)
4798 boundary_points[i][0] = 1. * (i + 1) / (n_points + 1);
4799
4800 const Quadrature<dim - 1> quadrature(boundary_points);
4801 const Quadrature<dim> q_projector(
4803 ReferenceCells::get_hypercube<dim>(), quadrature));
4804
4805 // next loop over all
4806 // boundary faces and
4807 // generate the info from
4808 // them
4809 for (const auto &cell : tria.active_cell_iterators())
4810 for (const unsigned int face_no :
4811 GeometryInfo<dim>::face_indices())
4812 {
4813 const typename ::Triangulation<dim, spacedim>::
4814 face_iterator face = cell->face(face_no);
4815
4816 if (face->at_boundary())
4817 {
4818 Point<dim> p0_dim(face->vertex(0));
4819 Point<2> p0(p0_dim[0], p0_dim[1]);
4820
4821 // loop over
4822 // all pieces
4823 // of the line
4824 // and generate
4825 // line-lets
4826 const auto offset =
4828 cell->reference_cell(),
4829 face_no,
4830 cell->combined_face_orientation(face_no),
4831 n_points);
4832 for (unsigned int i = 0; i < n_points; ++i)
4833 {
4834 const Point<dim> p1_dim(
4836 cell, q_projector.point(offset + i)));
4837 const Point<2> p1(p1_dim[0], p1_dim[1]);
4838
4839 line_list.emplace_back(p0,
4840 p1,
4841 face->user_flag_set(),
4842 cell->level());
4843 p0 = p1;
4844 }
4845
4846 // generate last piece
4847 const Point<dim> p1_dim(face->vertex(1));
4848 const Point<2> p1(p1_dim[0], p1_dim[1]);
4849 line_list.emplace_back(p0,
4850 p1,
4851 face->user_flag_set(),
4852 cell->level());
4853 }
4854 }
4855 }
4856
4857 break;
4858 }
4859
4860 case 3:
4861 {
4862 // curved boundary output
4863 // presently not supported
4864 Assert(mapping == nullptr, ExcNotImplemented());
4865
4866 // loop over all lines and compute their
4867 // projection on the plane perpendicular
4868 // to the direction of sight
4869
4870 // direction of view equals the unit
4871 // vector of the position of the
4872 // spectator to the origin.
4873 //
4874 // we chose here the viewpoint as in
4875 // gnuplot as default.
4876 //
4877 // TODO:[WB] Fix a potential problem with viewing angles in 3d Eps
4878 // GridOut
4879 // note: the following might be wrong
4880 // if one of the base vectors below
4881 // is in direction of the viewer, but
4882 // I am too tired at present to fix
4883 // this
4884 const double pi = numbers::PI;
4885 const double z_angle = eps_flags_3.azimut_angle;
4886 const double turn_angle = eps_flags_3.turn_angle;
4887 const Point<dim> view_direction(
4888 -std::sin(z_angle * 2. * pi / 360.) *
4889 std::sin(turn_angle * 2. * pi / 360.),
4890 +std::sin(z_angle * 2. * pi / 360.) *
4891 std::cos(turn_angle * 2. * pi / 360.),
4892 -std::cos(z_angle * 2. * pi / 360.));
4893
4894 // decide about the two unit vectors
4895 // in this plane. we chose the first one
4896 // to be the projection of the z-axis
4897 // to this plane
4898 const Tensor<1, dim> vector1 =
4899 Point<dim>(0, 0, 1) -
4900 ((Point<dim>(0, 0, 1) * view_direction) * view_direction);
4901 const Tensor<1, dim> unit_vector1 = vector1 / vector1.norm();
4902
4903 // now the third vector is fixed. we
4904 // chose the projection of a more or
4905 // less arbitrary vector to the plane
4906 // perpendicular to the first one
4907 const Tensor<1, dim> vector2 =
4908 (Point<dim>(1, 0, 0) -
4909 ((Point<dim>(1, 0, 0) * view_direction) * view_direction) -
4910 ((Point<dim>(1, 0, 0) * unit_vector1) * unit_vector1));
4911 const Tensor<1, dim> unit_vector2 = vector2 / vector2.norm();
4912
4913
4914 for (const auto &cell : tria.active_cell_iterators())
4915 for (const unsigned int line_no : cell->line_indices())
4916 {
4917 typename ::Triangulation<dim, spacedim>::line_iterator
4918 line = cell->line(line_no);
4919 line_list.emplace_back(
4920 Point<2>(line->vertex(0) * unit_vector2,
4921 line->vertex(0) * unit_vector1),
4922 Point<2>(line->vertex(1) * unit_vector2,
4923 line->vertex(1) * unit_vector1),
4924 line->user_flag_set(),
4925 cell->level());
4926 }
4927
4928 break;
4929 }
4930
4931 default:
4933 }
4934
4935
4936
4937 // find out minimum and maximum x and
4938 // y coordinates to compute offsets
4939 // and scaling factors
4940 double x_min = tria.begin_active()->vertex(0)[0];
4941 double x_max = x_min;
4942 double y_min = tria.begin_active()->vertex(0)[1];
4943 double y_max = y_min;
4944 unsigned int max_level = line_list.begin()->level;
4945
4946 for (LineList::const_iterator line = line_list.begin();
4947 line != line_list.end();
4948 ++line)
4949 {
4950 x_min = std::min(x_min, line->first[0]);
4951 x_min = std::min(x_min, line->second[0]);
4952
4953 x_max = std::max(x_max, line->first[0]);
4954 x_max = std::max(x_max, line->second[0]);
4955
4956 y_min = std::min(y_min, line->first[1]);
4957 y_min = std::min(y_min, line->second[1]);
4958
4959 y_max = std::max(y_max, line->first[1]);
4960 y_max = std::max(y_max, line->second[1]);
4961
4962 max_level = std::max(max_level, line->level);
4963 }
4964
4965 // scale in x-direction such that
4966 // in the output 0 <= x <= 300.
4967 // don't scale in y-direction to
4968 // preserve the shape of the
4969 // triangulation
4970 const double scale =
4971 (eps_flags_base.size /
4972 (eps_flags_base.size_type == GridOutFlags::EpsFlagsBase::width ?
4973 x_max - x_min :
4974 y_min - y_max));
4975
4976
4977 // now write preamble
4978 {
4979 // block this to have local
4980 // variables destroyed after
4981 // use
4982 std::time_t time1 = std::time(nullptr);
4983 std::tm *time = std::localtime(&time1);
4984 out << "%!PS-Adobe-2.0 EPSF-1.2" << '\n'
4985 << "%%Title: deal.II Output" << '\n'
4986 << "%%Creator: the deal.II library" << '\n'
4987 << "%%Creation Date: " << time->tm_year + 1900 << "/"
4988 << time->tm_mon + 1 << "/" << time->tm_mday << " - "
4989 << time->tm_hour << ":" << std::setw(2) << time->tm_min << ":"
4990 << std::setw(2) << time->tm_sec << '\n'
4991 << "%%BoundingBox: "
4992 // lower left corner
4993 << "0 0 "
4994 // upper right corner
4995 << static_cast<unsigned int>(
4996 std::floor(((x_max - x_min) * scale) + 1))
4997 << ' '
4998 << static_cast<unsigned int>(
4999 std::floor(((y_max - y_min) * scale) + 1))
5000 << '\n';
5001
5002 // define some abbreviations to keep
5003 // the output small:
5004 // m=move turtle to
5005 // x=execute line stroke
5006 // b=black pen
5007 // r=red pen
5008 out << "/m {moveto} bind def" << '\n'
5009 << "/x {lineto stroke} bind def" << '\n'
5010 << "/b {0 0 0 setrgbcolor} def" << '\n'
5011 << "/r {1 0 0 setrgbcolor} def" << '\n';
5012
5013 // calculate colors for level
5014 // coloring; level 0 is black,
5015 // other levels are blue
5016 // ... red
5017 if (eps_flags_base.color_lines_level)
5018 out << "/l { neg " << (max_level) << " add "
5019 << (0.66666 / std::max(1U, (max_level - 1)))
5020 << " mul 1 0.8 sethsbcolor} def" << '\n';
5021
5022 // in 2d, we can also plot cell
5023 // and vertex numbers, but this
5024 // requires a somewhat more
5025 // lengthy preamble. please
5026 // don't ask me what most of
5027 // this means, it is reverse
5028 // engineered from what GNUPLOT
5029 // uses in its output
5030 if ((dim == 2) && (eps_flags_2.write_cell_numbers ||
5031 eps_flags_2.write_vertex_numbers))
5032 {
5033 out
5034 << ("/R {rmoveto} bind def\n"
5035 "/Symbol-Oblique /Symbol findfont [1 0 .167 1 0 0] makefont\n"
5036 "dup length dict begin {1 index /FID eq {pop pop} {def} ifelse} forall\n"
5037 "currentdict end definefont\n"
5038 "/MFshow {{dup dup 0 get findfont exch 1 get scalefont setfont\n"
5039 "[ currentpoint ] exch dup 2 get 0 exch rmoveto dup dup 5 get exch 4 get\n"
5040 "{show} {stringwidth pop 0 rmoveto}ifelse dup 3 get\n"
5041 "{2 get neg 0 exch rmoveto pop} {pop aload pop moveto}ifelse} forall} bind def\n"
5042 "/MFwidth {0 exch {dup 3 get{dup dup 0 get findfont exch 1 get scalefont setfont\n"
5043 "5 get stringwidth pop add}\n"
5044 "{pop} ifelse} forall} bind def\n"
5045 "/MCshow { currentpoint stroke m\n"
5046 "exch dup MFwidth -2 div 3 -1 roll R MFshow } def\n")
5047 << '\n';
5048 }
5049
5050 out << "%%EndProlog" << '\n' << '\n';
5051
5052 // set fine lines
5053 out << eps_flags_base.line_width << " setlinewidth" << '\n';
5054 }
5055
5056 // now write the lines
5057 const Point<2> offset(x_min, y_min);
5058
5059 for (LineList::const_iterator line = line_list.begin();
5060 line != line_list.end();
5061 ++line)
5062 if (eps_flags_base.color_lines_level && (line->level > 0))
5063 out << line->level << " l " << (line->first - offset) * scale << " m "
5064 << (line->second - offset) * scale << " x" << '\n';
5065 else
5066 out << ((line->colorize && eps_flags_base.color_lines_on_user_flag) ?
5067 "r " :
5068 "b ")
5069 << (line->first - offset) * scale << " m "
5070 << (line->second - offset) * scale << " x" << '\n';
5071
5072 // finally write the cell numbers
5073 // in 2d, if that is desired
5074 if ((dim == 2) && (eps_flags_2.write_cell_numbers == true))
5075 {
5076 out << "(Helvetica) findfont 140 scalefont setfont" << '\n';
5077
5078 for (const auto &cell : tria.active_cell_iterators())
5079 {
5080 out << (cell->center()[0] - offset[0]) * scale << ' '
5081 << (cell->center()[1] - offset[1]) * scale << " m" << '\n'
5082 << "[ [(Helvetica) 12.0 0.0 true true (";
5083 if (eps_flags_2.write_cell_number_level)
5084 out << cell;
5085 else
5086 out << cell->index();
5087
5088 out << ")] "
5089 << "] -6 MCshow" << '\n';
5090 }
5091 }
5092
5093 // and the vertex numbers
5094 if ((dim == 2) && (eps_flags_2.write_vertex_numbers == true))
5095 {
5096 out << "(Helvetica) findfont 140 scalefont setfont" << '\n';
5097
5098 // have a list of those
5099 // vertices which we have
5100 // already tracked, to avoid
5101 // doing this multiply
5102 std::set<unsigned int> treated_vertices;
5103 for (const auto &cell : tria.active_cell_iterators())
5104 for (const unsigned int vertex_no : cell->vertex_indices())
5105 if (treated_vertices.find(cell->vertex_index(vertex_no)) ==
5106 treated_vertices.end())
5107 {
5108 treated_vertices.insert(cell->vertex_index(vertex_no));
5109
5110 out << (cell->vertex(vertex_no)[0] - offset[0]) * scale << ' '
5111 << (cell->vertex(vertex_no)[1] - offset[1]) * scale
5112 << " m" << '\n'
5113 << "[ [(Helvetica) 10.0 0.0 true true ("
5114 << cell->vertex_index(vertex_no) << ")] "
5115 << "] -6 MCshow" << '\n';
5116 }
5117 }
5118
5119 out << "showpage" << '\n';
5120
5121 // make sure everything now gets to
5122 // disk
5123 out.flush();
5124
5125 AssertThrow(out.fail() == false, ExcIO());
5126 }
5127 } // namespace
5128} // namespace internal
5129
5130
5131template <int dim, int spacedim>
5132void
5134 std::ostream &out,
5135 const Mapping<dim, spacedim> *mapping) const
5136{
5137 internal::write_eps(tria, out, mapping, eps_flags_2, eps_flags_3);
5138}
5139
5140
5141template <int dim, int spacedim>
5142void
5144 std::ostream &out,
5145 const OutputFormat output_format,
5146 const Mapping<dim, spacedim> *mapping) const
5147{
5148 switch (output_format)
5149 {
5150 case none:
5151 return;
5152
5153 case dx:
5154 write_dx(tria, out);
5155 return;
5156
5157 case ucd:
5158 write_ucd(tria, out);
5159 return;
5160
5161 case gnuplot:
5162 write_gnuplot(tria, out, mapping);
5163 return;
5164
5165 case eps:
5166 write_eps(tria, out, mapping);
5167 return;
5168
5169 case xfig:
5170 write_xfig(tria, out, mapping);
5171 return;
5172
5173 case msh:
5174 write_msh(tria, out);
5175 return;
5176
5177 case svg:
5178 write_svg(tria, out);
5179 return;
5180
5181 case mathgl:
5182 write_mathgl(tria, out);
5183 return;
5184
5185 case vtk:
5186 write_vtk(tria, out);
5187 return;
5188
5189 case vtu:
5190 write_vtu(tria, out);
5191 return;
5192 }
5193
5195}
5196
5197
5198template <int dim, int spacedim>
5199void
5201 std::ostream &out,
5202 const Mapping<dim, spacedim> *mapping) const
5203{
5204 write(tria, out, default_format, mapping);
5205}
5206
5207
5208// explicit instantiations
5209#include "grid/grid_out.inst"
5210
5211
void write_pvtu_record(std::ostream &out, const std::vector< std::string > &piece_names) const
void attach_triangulation(const Triangulation< dim, spacedim > &)
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation={})
virtual void build_patches(const unsigned int n_subdivisions=0)
Definition data_out.cc:1062
unsigned int n_boundary_faces(const Triangulation< dim, spacedim > &tria) const
Definition grid_out.cc:3786
GridOutFlags::Vtu vtu_flags
Definition grid_out.h:1615
GridOutFlags::Eps< 2 > eps_flags_2
Definition grid_out.h:1584
unsigned int write_ucd_lines(const Triangulation< dim, spacedim > &tria, const unsigned int next_element_index, std::ostream &out) const
Definition grid_out.cc:4107
void parse_parameters(ParameterHandler &param)
Definition grid_out.cc:700
void write_svg(const Triangulation< 2, 2 > &tria, std::ostream &out) const
Definition grid_out.cc:1702
void set_flags(const GridOutFlags::DX &flags)
Definition grid_out.cc:471
GridOutFlags::Eps< 1 > eps_flags_1
Definition grid_out.h:1578
GridOutFlags::XFig xfig_flags
Definition grid_out.h:1595
unsigned int write_ucd_faces(const Triangulation< dim, spacedim > &tria, const unsigned int next_element_index, std::ostream &out) const
Definition grid_out.cc:4065
void write_mathgl(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition grid_out.cc:2988
std::string default_suffix() const
Definition grid_out.cc:593
static void declare_parameters(ParameterHandler &param)
Definition grid_out.cc:651
static OutputFormat parse_output_format(const std::string &format_name)
Definition grid_out.cc:601
void write_vtk(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition grid_out.cc:3311
GridOutFlags::Gnuplot gnuplot_flags
Definition grid_out.h:1572
void write_msh(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition grid_out.cc:1021
unsigned int write_msh_lines(const Triangulation< dim, spacedim > &tria, const unsigned int next_element_index, std::ostream &out) const
Definition grid_out.cc:3945
void write_eps(const Triangulation< dim, spacedim > &tria, std::ostream &out, const Mapping< dim, spacedim > *mapping=nullptr) const
Definition grid_out.cc:5133
static std::string get_output_format_names()
Definition grid_out.cc:644
GridOutFlags::Eps< 3 > eps_flags_3
Definition grid_out.h:1590
void write(const Triangulation< dim, spacedim > &tria, std::ostream &out, const OutputFormat output_format, const Mapping< dim, spacedim > *mapping=nullptr) const
Definition grid_out.cc:5143
void write_vtu(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition grid_out.cc:3550
GridOutFlags::Ucd ucd_flags
Definition grid_out.h:1566
void write_dx(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition grid_out.cc:781
std::size_t memory_consumption() const
Definition grid_out.cc:746
unsigned int n_boundary_lines(const Triangulation< dim, spacedim > &tria) const
Definition grid_out.cc:3802
GridOutFlags::DX dx_flags
Definition grid_out.h:1554
unsigned int write_msh_faces(const Triangulation< dim, spacedim > &tria, const unsigned int next_element_index, std::ostream &out) const
Definition grid_out.cc:3906
void write_ucd(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition grid_out.cc:1126
GridOutFlags::Svg svg_flags
Definition grid_out.h:1600
OutputFormat default_format
Definition grid_out.h:1549
void write_xfig(const Triangulation< dim, spacedim > &tria, std::ostream &out, const Mapping< dim, spacedim > *mapping=nullptr) const
Definition grid_out.cc:1243
void write_gnuplot(const Triangulation< dim, spacedim > &tria, std::ostream &out, const Mapping< dim, spacedim > *mapping=nullptr) const
Definition grid_out.cc:4617
GridOutFlags::Vtk vtk_flags
Definition grid_out.h:1610
GridOutFlags::Msh msh_flags
Definition grid_out.h:1560
void write_mesh_per_processor_as_vtu(const Triangulation< dim, spacedim > &tria, const std::string &filename_without_extension, const bool view_levels=false, const bool include_artificial=false) const
Definition grid_out.cc:3597
@ vtk
write() calls write_vtk()
Definition grid_out.h:1021
@ eps
write() calls write_eps()
Definition grid_out.h:1009
@ msh
write() calls write_msh()
Definition grid_out.h:1015
@ xfig
write() calls write_xfig()
Definition grid_out.h:1013
@ dx
write() calls write_dx()
Definition grid_out.h:1005
@ ucd
write() calls write_ucd()
Definition grid_out.h:1011
@ gnuplot
write() calls write_gnuplot()
Definition grid_out.h:1007
@ mathgl
write() calls write_mathgl()
Definition grid_out.h:1019
@ svg
write() calls write_svg()
Definition grid_out.h:1017
@ none
Do nothing in write()
Definition grid_out.h:1003
@ vtu
write() calls write_vtu()
Definition grid_out.h:1023
GridOutFlags::MathGL mathgl_flags
Definition grid_out.h:1605
Abstract base class for mapping classes.
Definition mapping.h:320
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const =0
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const =0
void enter_subsection(const std::string &subsection, const bool create_path_if_needed=true)
long int get_integer(const std::string &entry_string) const
bool get_bool(const std::string &entry_name) const
void declare_entry(const std::string &entry, const std::string &default_value, const Patterns::PatternBase &pattern=Patterns::Anything(), const std::string &documentation="", const bool has_to_be_set=false)
std::string get(const std::string &entry_string) const
double get_double(const std::string &entry_name) const
Definition point.h:113
static DataSetDescriptor face(const ReferenceCell &reference_cell, const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points)
Class which transforms dim - 1-dimensional quadrature rules to dim-dimensional face quadratures.
Definition qprojector.h:69
static Quadrature< dim > project_to_all_faces(const ReferenceCell &reference_cell, const hp::QCollection< dim - 1 > &quadrature)
const Point< dim > & point(const unsigned int i) const
unsigned int gmsh_element_type() const
numbers::NumberTraits< Number >::real_type norm() const
void save_user_flags_line(std::ostream &out) const
void save(Archive &ar, const unsigned int version) const
cell_iterator begin(const unsigned int level=0) const
virtual types::subdomain_id locally_owned_subdomain() const
unsigned int n_active_cells() const
const std::vector< Point< spacedim > > & get_vertices() const
unsigned int n_used_vertices() const
const std::vector< ReferenceCell > & get_reference_cells() const
cell_iterator end() const
const std::vector< bool > & get_used_vertices() const
active_cell_iterator begin_active(const unsigned int level=0) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
Point< 2 > second
Definition grid_out.cc:4633
bool colorize
Definition grid_out.cc:4634
Point< 2 > first
Definition grid_out.cc:4632
unsigned int level
Definition grid_out.cc:4635
unsigned int vertex_indices[2]
const unsigned int v0
const unsigned int v1
IteratorRange< active_face_iterator > active_face_iterators() const
IteratorRange< active_cell_iterator > active_cell_iterators() const
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
IteratorRange< cell_iterator > cell_iterators() const
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcInvalidState()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
std::size_t size
Definition mpi.cc:745
const unsigned int n_procs
Definition mpi.cc:935
void write_eps(const std::vector< Patch< 2, spacedim > > &patches, const std::vector< std::string > &data_names, const std::vector< std::tuple< unsigned int, unsigned int, std::string, DataComponentInterpretation::DataComponentInterpretation > > &nonscalar_data_ranges, const EpsFlags &flags, std::ostream &out)
void write_vtu_header(std::ostream &out, const VtkFlags &flags)
void write_vtu(const std::vector< Patch< dim, spacedim > > &patches, const std::vector< std::string > &data_names, const std::vector< std::tuple< unsigned int, unsigned int, std::string, DataComponentInterpretation::DataComponentInterpretation > > &nonscalar_data_ranges, const VtkFlags &flags, std::ostream &out)
void write_vtu_main(const std::vector< Patch< dim, spacedim > > &patches, const std::vector< std::string > &data_names, const std::vector< std::tuple< unsigned int, unsigned int, std::string, DataComponentInterpretation::DataComponentInterpretation > > &nonscalar_data_ranges, const VtkFlags &flags, std::ostream &out)
void write_vtu_footer(std::ostream &out)
void write_gnuplot(const std::vector< Patch< dim, spacedim > > &patches, const std::vector< std::string > &data_names, const std::vector< std::tuple< unsigned int, unsigned int, std::string, DataComponentInterpretation::DataComponentInterpretation > > &nonscalar_data_ranges, const GnuplotFlags &flags, std::ostream &out)
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:193
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr ReferenceCell Triangle
constexpr ReferenceCell Hexahedron
constexpr ReferenceCell Pyramid
constexpr ReferenceCell Wedge
constexpr ReferenceCell Vertex
constexpr ReferenceCell Quadrilateral
constexpr ReferenceCell Tetrahedron
constexpr ReferenceCell Line
VectorType::value_type * end(VectorType &V)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:99
std::string encode_base64(const std::vector< unsigned char > &binary_input)
Definition utilities.cc:427
std::string compress(const std::string &input)
Definition utilities.cc:383
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:466
constexpr double PI
Definition numbers.h:240
constexpr unsigned int invalid_unsigned_int
Definition types.h:238
constexpr types::boundary_id internal_face_boundary_id
Definition types.h:329
constexpr types::boundary_id invalid_boundary_id
Definition types.h:309
constexpr types::manifold_id flat_manifold_id
Definition types.h:342
constexpr types::subdomain_id artificial_subdomain_id
Definition types.h:402
typename type_identity< T >::type type_identity_t
Definition type_traits.h:95
::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:184
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