Reference documentation for deal.II version GIT relicensing-136-gb80d0be4af 2024-03-18 08:20:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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") == std::string("width"))
235 else if (param.get("Size by") == std::string("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 viw point, that is, the angle "
348 "in the plane from the x-axis.");
349 param.declare_entry("Elevation",
350 "30",
352 "Elevation of the view point above the xy-plane.");
353 }
354
355
356 void
358 {
360 azimut_angle = 90 - param.get_double("Elevation");
361 turn_angle = param.get_double("Azimuth");
362 }
363
364
365
367 : draw_boundary(true)
368 , color_by(material_id)
369 , level_depth(true)
370 , n_boundary_face_points(0)
371 , scaling(1., 1.)
372 , fill_style(20)
373 , line_style(0)
374 , line_thickness(1)
375 , boundary_style(0)
376 , boundary_thickness(3)
377 {}
378
379
380 void
382 {
383 param.declare_entry("Boundary", "true", Patterns::Bool());
384 param.declare_entry("Level color", "false", Patterns::Bool());
385 param.declare_entry("Level depth", "true", Patterns::Bool());
386 // TODO: Unify this number with other output formats
387 param.declare_entry("Boundary points", "0", Patterns::Integer());
388 param.declare_entry("Fill style", "20", Patterns::Integer());
389 param.declare_entry("Line style", "0", Patterns::Integer());
390 param.declare_entry("Line width", "1", Patterns::Integer());
391 param.declare_entry("Boundary style", "0", Patterns::Integer());
392 param.declare_entry("Boundary width", "3", Patterns::Integer());
393 }
394
395
396 void
398 {
399 draw_boundary = param.get_bool("Boundary");
400 level_depth = param.get_bool("Level depth");
401 n_boundary_face_points = param.get_integer("Boundary points");
402 fill_style = param.get_integer("Fill style");
403 line_style = param.get_integer("Line style");
404 line_thickness = param.get_integer("Line width");
405 boundary_style = param.get_integer("Boundary style");
406 boundary_thickness = param.get_integer("Boundary width");
407 }
408
409 Svg::Svg(const unsigned int line_thickness,
410 const unsigned int boundary_line_thickness,
411 bool margin,
412 const Background background,
413 const int azimuth_angle,
414 const int polar_angle,
415 const Coloring coloring,
416 const bool convert_level_number_to_height,
417 const bool label_level_number,
418 const bool label_cell_index,
419 const bool label_material_id,
420 const bool label_subdomain_id,
421 const bool draw_colorbar,
422 const bool draw_legend,
423 const bool label_boundary_id)
424 : height(1000)
425 , width(0)
426 , line_thickness(line_thickness)
427 , boundary_line_thickness(boundary_line_thickness)
428 , margin(margin)
429 , background(background)
430 , azimuth_angle(azimuth_angle)
431 , polar_angle(polar_angle)
432 , coloring(coloring)
433 , convert_level_number_to_height(convert_level_number_to_height)
434 , level_height_factor(0.3f)
435 , cell_font_scaling(1.f)
436 , label_level_number(label_level_number)
437 , label_cell_index(label_cell_index)
438 , label_material_id(label_material_id)
439 , label_subdomain_id(label_subdomain_id)
440 , label_level_subdomain_id(false)
441 , label_boundary_id(label_boundary_id)
442 , draw_colorbar(draw_colorbar)
443 , draw_legend(draw_legend)
444 {}
445
447 : draw_bounding_box(false) // box
448 {}
449
450 void
452 {
453 param.declare_entry("Draw bounding box", "false", Patterns::Bool());
454 }
455
456 void
458 {
459 draw_bounding_box = param.get_bool("Draw bounding box");
460 }
461} // end namespace GridOutFlags
462
463
464
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 h = .6 - (index / (n - 1.)) * .6;
2143 else
2144 h = .6;
2145
2146 unsigned int r = 0;
2147 unsigned int g = 0;
2148 unsigned int b = 0;
2149
2150 unsigned int i = static_cast<unsigned int>(h * 6);
2151
2152 double f = h * 6 - i;
2153 double q = 1 - f;
2154 double t = f;
2155
2156 switch (i % 6)
2157 {
2158 case 0:
2159 r = 255, g = static_cast<unsigned int>(.5 + 255 * t);
2160 break;
2161 case 1:
2162 r = static_cast<unsigned int>(.5 + 255 * q), g = 255;
2163 break;
2164 case 2:
2165 g = 255, b = static_cast<unsigned int>(.5 + 255 * t);
2166 break;
2167 case 3:
2168 g = static_cast<unsigned int>(.5 + 255 * q), b = 255;
2169 break;
2170 case 4:
2171 r = static_cast<unsigned int>(.5 + 255 * t), b = 255;
2172 break;
2173 case 5:
2174 r = 255, b = static_cast<unsigned int>(.5 + 255 * q);
2175 break;
2176 default:
2177 break;
2178 }
2179
2180 switch (svg_flags.coloring)
2181 {
2183 labeling_index = *materials_it++;
2184 break;
2186 labeling_index = *levels_it++;
2187 break;
2189 labeling_index = *subdomains_it++;
2190 break;
2192 labeling_index = *level_subdomains_it++;
2193 break;
2194 default:
2195 break;
2196 }
2197
2198 out << " path.p" << labeling_index << "{fill:rgb(" << r << ',' << g
2199 << ',' << b << "); "
2200 << "stroke:rgb(25,25,25); stroke-width:"
2201 << svg_flags.line_thickness << '}' << '\n';
2202
2203 out << " path.ps" << labeling_index << "{fill:rgb("
2204 << static_cast<unsigned int>(.5 + .75 * r) << ','
2205 << static_cast<unsigned int>(.5 + .75 * g) << ','
2206 << static_cast<unsigned int>(.5 + .75 * b) << "); "
2207 << "stroke:rgb(20,20,20); stroke-width:"
2208 << svg_flags.line_thickness << '}' << '\n';
2209
2210 out << " rect.r" << labeling_index << "{fill:rgb(" << r << ',' << g
2211 << ',' << b << "); "
2212 << "stroke:rgb(25,25,25); stroke-width:"
2213 << svg_flags.line_thickness << '}' << '\n';
2214
2215 ++labeling_index;
2216 }
2217 }
2218
2219 out << "]]></style>" << '\n' << '\n';
2220
2221 // background rectangle
2222 out << " <rect class=\"background\" width=\"" << width << "\" height=\""
2223 << height << "\"/>" << '\n';
2224
2226 {
2227 unsigned int x_offset = 0;
2228
2229 if (svg_flags.margin)
2230 x_offset = static_cast<unsigned int>(.5 + (height / 100.) *
2231 (margin_in_percent / 2.));
2232 else
2233 x_offset = static_cast<unsigned int>(.5 + height * .025);
2234
2235 out
2236 << " <text x=\"" << x_offset << "\" y=\""
2237 << static_cast<unsigned int>(.5 + height * .0525) << '\"'
2238 << " style=\"font-weight:100; fill:lightsteelblue; text-anchor:start; font-family:Courier; font-size:"
2239 << static_cast<unsigned int>(.5 + height * .045) << "px\">"
2240 << "deal.II"
2241 << "</text>" << '\n';
2242
2243 // out << " <text x=\"" << x_offset + static_cast<unsigned int>(.5 +
2244 // height * .045 * 4.75) << "\" y=\"" << static_cast<unsigned int>(.5 +
2245 // height * .0525) << '\"'
2246 // << " style=\"fill:lightsteelblue; text-anchor:start; font-size:" <<
2247 // font_size << "\">"
2248 // << now->tm_mday << '/' << now->tm_mon + 1 << '/' << now->tm_year +
2249 // 1900
2250 // << " - " << now->tm_hour << ':';
2251 //
2252 // if(now->tm_min < 10) out << '0';
2253 //
2254 // out << now->tm_min
2255 // << "</text>"<< '\n' << '\n';
2256 }
2257
2258 // draw the cells, starting out from the minimal level (in order to guaranty a
2259 // correct perspective view)
2260 out << " <!-- cells -->" << '\n';
2261
2262 for (unsigned int level_index = min_level; level_index <= max_level;
2263 level_index++)
2264 {
2265 for (const auto &cell : tria.cell_iterators_on_level(level_index))
2266 {
2267 if (!svg_flags.convert_level_number_to_height && !cell->is_active())
2268 continue;
2269
2270 // draw the current cell
2271 out << " <path";
2272
2273 if (svg_flags.coloring != 0u)
2274 {
2275 out << " class=\"p";
2276
2277 if (!cell->is_active() &&
2279 out << 's';
2280
2281 switch (svg_flags.coloring)
2282 {
2284 out << cell->material_id();
2285 break;
2287 out << static_cast<unsigned int>(cell->level());
2288 break;
2290 if (cell->is_active())
2291 out << cell->subdomain_id() + 2;
2292 else
2293 out << 'X';
2294 break;
2296 out << cell->level_subdomain_id() + 2;
2297 break;
2298 default:
2299 break;
2300 }
2301
2302 out << '\"';
2303 }
2304
2305 out << " d=\"M ";
2306
2307 point[0] = cell->vertex(0)[0];
2308 point[1] = cell->vertex(0)[1];
2309 point[2] = 0;
2310
2312 {
2313 point[2] = svg_flags.level_height_factor *
2314 (static_cast<float>(cell->level()) /
2315 static_cast<float>(n_levels)) *
2316 std::max(x_dimension, y_dimension);
2317 }
2318
2319 projection_decomposition = svg_project_point(point,
2320 camera_position,
2321 camera_direction,
2322 camera_horizontal,
2323 camera_focus);
2324
2325 out << static_cast<unsigned int>(
2326 .5 +
2327 ((projection_decomposition[0] - x_min_perspective) /
2328 x_dimension_perspective) *
2329 (width - (width / 100.) * 2. * margin_in_percent) +
2330 ((width / 100.) * margin_in_percent))
2331 << ' '
2332 << static_cast<unsigned int>(
2333 .5 + height - (height / 100.) * margin_in_percent -
2334 ((projection_decomposition[1] - y_min_perspective) /
2335 y_dimension_perspective) *
2336 (height - (height / 100.) * 2. * margin_in_percent));
2337
2338 out << " L ";
2339
2340 point[0] = cell->vertex(1)[0];
2341 point[1] = cell->vertex(1)[1];
2342
2343 projection_decomposition = svg_project_point(point,
2344 camera_position,
2345 camera_direction,
2346 camera_horizontal,
2347 camera_focus);
2348
2349 out << static_cast<unsigned int>(
2350 .5 +
2351 ((projection_decomposition[0] - x_min_perspective) /
2352 x_dimension_perspective) *
2353 (width - (width / 100.) * 2. * margin_in_percent) +
2354 ((width / 100.) * margin_in_percent))
2355 << ' '
2356 << static_cast<unsigned int>(
2357 .5 + height - (height / 100.) * margin_in_percent -
2358 ((projection_decomposition[1] - y_min_perspective) /
2359 y_dimension_perspective) *
2360 (height - (height / 100.) * 2. * margin_in_percent));
2361
2362 out << " L ";
2363
2364 if (cell->n_vertices() == 4) // in case of quadrilateral
2365 {
2366 point[0] = cell->vertex(3)[0];
2367 point[1] = cell->vertex(3)[1];
2368
2369 projection_decomposition = svg_project_point(point,
2370 camera_position,
2371 camera_direction,
2372 camera_horizontal,
2373 camera_focus);
2374
2375 out << static_cast<unsigned int>(
2376 .5 +
2377 ((projection_decomposition[0] - x_min_perspective) /
2378 x_dimension_perspective) *
2379 (width - (width / 100.) * 2. * margin_in_percent) +
2380 ((width / 100.) * margin_in_percent))
2381 << ' '
2382 << static_cast<unsigned int>(
2383 .5 + height - (height / 100.) * margin_in_percent -
2384 ((projection_decomposition[1] - y_min_perspective) /
2385 y_dimension_perspective) *
2386 (height - (height / 100.) * 2. * margin_in_percent));
2387
2388 out << " L ";
2389 }
2390
2391 point[0] = cell->vertex(2)[0];
2392 point[1] = cell->vertex(2)[1];
2393
2394 projection_decomposition = svg_project_point(point,
2395 camera_position,
2396 camera_direction,
2397 camera_horizontal,
2398 camera_focus);
2399
2400 out << static_cast<unsigned int>(
2401 .5 +
2402 ((projection_decomposition[0] - x_min_perspective) /
2403 x_dimension_perspective) *
2404 (width - (width / 100.) * 2. * margin_in_percent) +
2405 ((width / 100.) * margin_in_percent))
2406 << ' '
2407 << static_cast<unsigned int>(
2408 .5 + height - (height / 100.) * margin_in_percent -
2409 ((projection_decomposition[1] - y_min_perspective) /
2410 y_dimension_perspective) *
2411 (height - (height / 100.) * 2. * margin_in_percent));
2412
2413 out << " L ";
2414
2415 point[0] = cell->vertex(0)[0];
2416 point[1] = cell->vertex(0)[1];
2417
2418 projection_decomposition = svg_project_point(point,
2419 camera_position,
2420 camera_direction,
2421 camera_horizontal,
2422 camera_focus);
2423
2424 out << static_cast<unsigned int>(
2425 .5 +
2426 ((projection_decomposition[0] - x_min_perspective) /
2427 x_dimension_perspective) *
2428 (width - (width / 100.) * 2. * margin_in_percent) +
2429 ((width / 100.) * margin_in_percent))
2430 << ' '
2431 << static_cast<unsigned int>(
2432 .5 + height - (height / 100.) * margin_in_percent -
2433 ((projection_decomposition[1] - y_min_perspective) /
2434 y_dimension_perspective) *
2435 (height - (height / 100.) * 2. * margin_in_percent));
2436
2437 out << "\"/>" << '\n';
2438
2439 // label the current cell
2443 {
2444 point[0] = cell->center()[0];
2445 point[1] = cell->center()[1];
2446 point[2] = 0;
2447
2449 {
2450 point[2] = svg_flags.level_height_factor *
2451 (static_cast<float>(cell->level()) /
2452 static_cast<float>(n_levels)) *
2453 std::max(x_dimension, y_dimension);
2454 }
2455
2456 const double distance_to_camera =
2457 std::hypot(point[0] - camera_position[0],
2458 point[1] - camera_position[1],
2459 point[2] - camera_position[2]);
2460 const double distance_factor =
2461 distance_to_camera / (2. * std::max(x_dimension, y_dimension));
2462
2463 projection_decomposition = svg_project_point(point,
2464 camera_position,
2465 camera_direction,
2466 camera_horizontal,
2467 camera_focus);
2468
2469 const unsigned int font_size_this_cell =
2470 static_cast<unsigned int>(
2471 .5 +
2472 cell_label_font_size *
2473 std::pow(.5, cell->level() - 4. + 3.5 * distance_factor));
2474
2475 out << " <text"
2476 << " x=\""
2477 << static_cast<unsigned int>(
2478 .5 +
2479 ((projection_decomposition[0] - x_min_perspective) /
2480 x_dimension_perspective) *
2481 (width - (width / 100.) * 2. * margin_in_percent) +
2482 ((width / 100.) * margin_in_percent))
2483 << "\" y=\""
2484 << static_cast<unsigned int>(
2485 .5 + height - (height / 100.) * margin_in_percent -
2486 ((projection_decomposition[1] - y_min_perspective) /
2487 y_dimension_perspective) *
2488 (height - (height / 100.) * 2. * margin_in_percent) +
2489 0.5 * font_size_this_cell)
2490 << "\" style=\"font-size:" << font_size_this_cell << "px\">";
2491
2493 {
2494 out << cell->level();
2495 }
2496
2498 {
2500 out << '.';
2501 out << cell->index();
2502 }
2503
2505 {
2508 out << ',';
2509 out << static_cast<std::make_signed_t<types::material_id>>(
2510 cell->material_id());
2511 }
2512
2514 {
2517 out << ',';
2518 if (cell->is_active())
2519 out << static_cast<std::make_signed_t<types::subdomain_id>>(
2520 cell->subdomain_id());
2521 else
2522 out << 'X';
2523 }
2524
2526 {
2531 out << ',';
2532 out << static_cast<std::make_signed_t<types::subdomain_id>>(
2533 cell->level_subdomain_id());
2534 }
2535
2536 out << "</text>" << '\n';
2537 }
2538
2539 // if the current cell lies at the boundary of the triangulation, draw
2540 // the additional boundary line
2542 {
2543 for (auto faceIndex : cell->face_indices())
2544 {
2545 if (cell->at_boundary(faceIndex))
2546 {
2547 point[0] = cell->face(faceIndex)->vertex(0)[0];
2548 point[1] = cell->face(faceIndex)->vertex(0)[1];
2549 point[2] = 0;
2550
2552 {
2553 point[2] = svg_flags.level_height_factor *
2554 (static_cast<float>(cell->level()) /
2555 static_cast<float>(n_levels)) *
2556 std::max(x_dimension, y_dimension);
2557 }
2558
2559 projection_decomposition =
2560 svg_project_point(point,
2561 camera_position,
2562 camera_direction,
2563 camera_horizontal,
2564 camera_focus);
2565
2566 out << " <line x1=\""
2567 << static_cast<unsigned int>(
2568 .5 +
2569 ((projection_decomposition[0] -
2570 x_min_perspective) /
2571 x_dimension_perspective) *
2572 (width -
2573 (width / 100.) * 2. * margin_in_percent) +
2574 ((width / 100.) * margin_in_percent))
2575 << "\" y1=\""
2576 << static_cast<unsigned int>(
2577 .5 + height -
2578 (height / 100.) * margin_in_percent -
2579 ((projection_decomposition[1] -
2580 y_min_perspective) /
2581 y_dimension_perspective) *
2582 (height -
2583 (height / 100.) * 2. * margin_in_percent));
2584
2585 point[0] = cell->face(faceIndex)->vertex(1)[0];
2586 point[1] = cell->face(faceIndex)->vertex(1)[1];
2587 point[2] = 0;
2588
2590 {
2591 point[2] = svg_flags.level_height_factor *
2592 (static_cast<float>(cell->level()) /
2593 static_cast<float>(n_levels)) *
2594 std::max(x_dimension, y_dimension);
2595 }
2596
2597 projection_decomposition =
2598 svg_project_point(point,
2599 camera_position,
2600 camera_direction,
2601 camera_horizontal,
2602 camera_focus);
2603
2604 out << "\" x2=\""
2605 << static_cast<unsigned int>(
2606 .5 +
2607 ((projection_decomposition[0] -
2608 x_min_perspective) /
2609 x_dimension_perspective) *
2610 (width -
2611 (width / 100.) * 2. * margin_in_percent) +
2612 ((width / 100.) * margin_in_percent))
2613 << "\" y2=\""
2614 << static_cast<unsigned int>(
2615 .5 + height -
2616 (height / 100.) * margin_in_percent -
2617 ((projection_decomposition[1] -
2618 y_min_perspective) /
2619 y_dimension_perspective) *
2620 (height -
2621 (height / 100.) * 2. * margin_in_percent))
2622 << "\"/>" << '\n';
2623
2624
2626 {
2627 const double distance_to_camera =
2628 std::hypot(point[0] - camera_position[0],
2629 point[1] - camera_position[1],
2630 point[2] - camera_position[2]);
2631 const double distance_factor =
2632 distance_to_camera /
2633 (2. * std::max(x_dimension, y_dimension));
2634
2635 const unsigned int font_size_this_edge =
2636 static_cast<unsigned int>(
2637 .5 + .5 * cell_label_font_size *
2638 std::pow(.5,
2639 cell->level() - 4. +
2640 3.5 * distance_factor));
2641
2642 point[0] = cell->face(faceIndex)->center()[0];
2643 point[1] = cell->face(faceIndex)->center()[1];
2644 point[2] = 0;
2645
2647 {
2648 point[2] = svg_flags.level_height_factor *
2649 (static_cast<float>(cell->level()) /
2650 static_cast<float>(n_levels)) *
2651 std::max(x_dimension, y_dimension);
2652 }
2653
2654 projection_decomposition =
2655 svg_project_point(point,
2656 camera_position,
2657 camera_direction,
2658 camera_horizontal,
2659 camera_focus);
2660
2661 const unsigned int xc = static_cast<unsigned int>(
2662 .5 +
2663 ((projection_decomposition[0] - x_min_perspective) /
2664 x_dimension_perspective) *
2665 (width -
2666 (width / 100.) * 2. * margin_in_percent) +
2667 ((width / 100.) * margin_in_percent));
2668 const unsigned int yc = static_cast<unsigned int>(
2669 .5 + height - (height / 100.) * margin_in_percent -
2670 ((projection_decomposition[1] - y_min_perspective) /
2671 y_dimension_perspective) *
2672 (height -
2673 (height / 100.) * 2. * margin_in_percent));
2674
2675 out << " <circle cx=\"" << xc << "\" cy=\"" << yc
2676 << "\" r=\"" << font_size_this_edge << "\" />"
2677 << '\n';
2678
2679 out << " <text x=\"" << xc << "\" y=\"" << yc
2680 << "\" style=\"font-size:" << font_size_this_edge
2681 << "px\" dominant-baseline=\"middle\">"
2682 << static_cast<int>(
2683 cell->face(faceIndex)->boundary_id())
2684 << "</text>" << '\n';
2685 }
2686 }
2687 }
2688 }
2689 }
2690 }
2691
2692
2693
2694 // draw the legend
2696 out << '\n' << " <!-- legend -->" << '\n';
2697
2698 additional_width = 0;
2699 if (!svg_flags.margin)
2700 additional_width = static_cast<unsigned int>(.5 + (height / 100.) * 2.5);
2701
2702 // explanation of the cell labeling
2703 if (svg_flags.draw_legend &&
2707 {
2708 unsigned int line_offset = 0;
2709 out << " <rect x=\"" << width + additional_width << "\" y=\""
2710 << static_cast<unsigned int>(.5 + (height / 100.) * margin_in_percent)
2711 << "\" width=\""
2712 << static_cast<unsigned int>(.5 + (height / 100.) *
2713 (40. - margin_in_percent))
2714 << "\" height=\"" << static_cast<unsigned int>(.5 + height * .215)
2715 << "\"/>" << '\n';
2716
2717 out << " <text x=\""
2718 << width + additional_width +
2719 static_cast<unsigned int>(.5 + (height / 100.) * 1.25)
2720 << "\" y=\""
2721 << static_cast<unsigned int>(.5 +
2722 (height / 100.) * margin_in_percent +
2723 (++line_offset) * 1.5 * font_size)
2724 << "\" style=\"text-anchor:start; font-weight:bold; font-size:"
2725 << font_size << "px\">"
2726 << "cell label"
2727 << "</text>" << '\n';
2728
2730 {
2731 out << " <text x=\""
2732 << width + additional_width +
2733 static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2734 << "\" y=\""
2735 << static_cast<unsigned int>(.5 +
2736 (height / 100.) * margin_in_percent +
2737 (++line_offset) * 1.5 * font_size)
2738 << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2739 << font_size << "px\">"
2740 << "cell_level";
2741
2745 out << '.';
2746
2747 out << "</text>" << '\n';
2748 }
2749
2751 {
2752 out << " <text x=\""
2753 << width + additional_width +
2754 static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2755 << "\" y=\""
2756 << static_cast<unsigned int>(.5 +
2757 (height / 100.) * margin_in_percent +
2758 (++line_offset) * 1.5 * font_size)
2759 << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2760 << font_size << "px\">"
2761 << "cell_index";
2762
2765 out << ',';
2766
2767 out << "</text>" << '\n';
2768 }
2769
2771 {
2772 out << " <text x=\""
2773 << width + additional_width +
2774 static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2775 << "\" y=\""
2776 << static_cast<unsigned int>(.5 +
2777 (height / 100.) * margin_in_percent +
2778 (++line_offset) * 1.5 * font_size)
2779 << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2780 << font_size << "px\">"
2781 << "material_id";
2782
2785 out << ',';
2786
2787 out << "</text>" << '\n';
2788 }
2789
2791 {
2792 out << " <text x= \""
2793 << width + additional_width +
2794 static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2795 << "\" y=\""
2796 << static_cast<unsigned int>(.5 +
2797 (height / 100.) * margin_in_percent +
2798 (++line_offset) * 1.5 * font_size)
2799 << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2800 << font_size << "px\">"
2801 << "subdomain_id";
2802
2804 out << ',';
2805
2806 out << "</text>" << '\n';
2807 }
2808
2810 {
2811 out << " <text x= \""
2812 << width + additional_width +
2813 static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2814 << "\" y=\""
2815 << static_cast<unsigned int>(.5 +
2816 (height / 100.) * margin_in_percent +
2817 (++line_offset) * 1.5 * font_size)
2818 << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2819 << font_size << "px\">"
2820 << "level_subdomain_id"
2821 << "</text>" << '\n';
2822 }
2823
2825 {
2826 out << " <text x=\""
2827 << width + additional_width +
2828 static_cast<unsigned int>(.5 + (height / 100.) * 1.25)
2829 << "\" y=\""
2830 << static_cast<unsigned int>(.5 +
2831 (height / 100.) * margin_in_percent +
2832 (++line_offset) * 1.5 * font_size)
2833 << "\" style=\"text-anchor:start; font-weight:bold; font-size:"
2834 << font_size << "px\">"
2835 << "edge label"
2836 << "</text>" << '\n';
2837
2838 out << " <text x= \""
2839 << width + additional_width +
2840 static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2841 << "\" y=\""
2842 << static_cast<unsigned int>(.5 +
2843 (height / 100.) * margin_in_percent +
2844 (++line_offset) * 1.5 * font_size)
2845 << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2846 << font_size << "px\">"
2847 << "boundary_id"
2848 << "</text>" << '\n';
2849 }
2850 }
2851
2852 // show azimuth angle and polar angle as text below the explanation of the
2853 // cell labeling
2855 {
2856 out << " <text x=\"" << width + additional_width << "\" y=\""
2857 << static_cast<unsigned int>(
2858 .5 + (height / 100.) * margin_in_percent + 13.75 * font_size)
2859 << "\" style=\"text-anchor:start; font-size:" << font_size << "px\">"
2860 << "azimuth: " << svg_flags.azimuth_angle
2861 << "°, polar: " << svg_flags.polar_angle << "°</text>" << '\n';
2862 }
2863
2864
2865 // draw the colorbar
2867 {
2868 out << '\n' << " <!-- colorbar -->" << '\n';
2869
2870 out << " <text x=\"" << width + additional_width << "\" y=\""
2871 << static_cast<unsigned int>(
2872 .5 + (height / 100.) * (margin_in_percent + 29.) -
2873 (font_size / 1.25))
2874 << "\" style=\"text-anchor:start; font-weight:bold; font-size:"
2875 << font_size << "px\">";
2876
2877 switch (svg_flags.coloring)
2878 {
2879 case 1:
2880 out << "material_id";
2881 break;
2882 case 2:
2883 out << "level_number";
2884 break;
2885 case 3:
2886 out << "subdomain_id";
2887 break;
2888 case 4:
2889 out << "level_subdomain_id";
2890 break;
2891 default:
2892 break;
2893 }
2894
2895 out << "</text>" << '\n';
2896
2897 unsigned int element_height = static_cast<unsigned int>(
2898 ((height / 100.) * (71. - 2. * margin_in_percent)) / n);
2899 unsigned int element_width =
2900 static_cast<unsigned int>(.5 + (height / 100.) * 2.5);
2901
2902 int labeling_index = 0;
2903 auto materials_it = materials.begin();
2904 auto levels_it = levels.begin();
2905 auto subdomains_it = subdomains.begin();
2906 auto level_subdomains_it = level_subdomains.begin();
2907
2908 for (unsigned int index = 0; index < n; ++index)
2909 {
2910 switch (svg_flags.coloring)
2911 {
2913 labeling_index = *materials_it++;
2914 break;
2916 labeling_index = *levels_it++;
2917 break;
2919 labeling_index = *subdomains_it++;
2920 break;
2922 labeling_index = *level_subdomains_it++;
2923 break;
2924 default:
2925 break;
2926 }
2927
2928 out << " <rect class=\"r" << labeling_index << "\" x=\""
2929 << width + additional_width << "\" y=\""
2930 << static_cast<unsigned int>(.5 + (height / 100.) *
2931 (margin_in_percent + 29)) +
2932 (n - index - 1) * element_height
2933 << "\" width=\"" << element_width << "\" height=\""
2934 << element_height << "\"/>" << '\n';
2935
2936 out << " <text x=\""
2937 << width + additional_width + 1.5 * element_width << "\" y=\""
2938 << static_cast<unsigned int>(.5 + (height / 100.) *
2939 (margin_in_percent + 29)) +
2940 (n - index - 1 + .5) * element_height +
2941 static_cast<unsigned int>(.5 + font_size * .35)
2942 << "\""
2943 << " style=\"text-anchor:start; font-size:"
2944 << static_cast<unsigned int>(.5 + font_size) << "px";
2945
2946 if (index == 0 || index == n - 1)
2947 out << "; font-weight:bold";
2948
2949 out << "\">" << labeling_index;
2950
2951 if (index == n - 1)
2952 out << " max";
2953 if (index == 0)
2954 out << " min";
2955
2956 out << "</text>" << '\n';
2957
2958 ++labeling_index;
2959 }
2960 }
2961
2962
2963 // finalize the svg file
2964 out << '\n' << "</svg>";
2965 out.flush();
2966}
2967
2968
2969
2970template <>
2971void
2972GridOut::write_mathgl(const Triangulation<1> &, std::ostream &) const
2973{
2974 // 1d specialization not done yet
2976}
2977
2978
2979
2980template <int dim, int spacedim>
2981void
2983 std::ostream &out) const
2984{
2985 AssertThrow(out.fail() == false, ExcIO());
2986
2987 // (i) write header
2988 {
2989 // block this to have local variables destroyed after use
2990 const std::time_t time1 = std::time(nullptr);
2991 const std::tm *time = std::localtime(&time1);
2992
2993 out
2994 << "\n#"
2995 << "\n# This file was generated by the deal.II library."
2996 << "\n# Date = " << time->tm_year + 1900 << "/" << std::setfill('0')
2997 << std::setw(2) << time->tm_mon + 1 << "/" << std::setfill('0')
2998 << std::setw(2) << time->tm_mday << "\n# Time = " << std::setfill('0')
2999 << std::setw(2) << time->tm_hour << ":" << std::setfill('0')
3000 << std::setw(2) << time->tm_min << ":" << std::setfill('0')
3001 << std::setw(2) << time->tm_sec << "\n#"
3002 << "\n# For a description of the MathGL script format see the MathGL manual. "
3003 << "\n#"
3004 << "\n# Note: This file is understood by MathGL v2.1 and higher only, and can "
3005 << "\n# be quickly viewed in a graphical environment using \'mglview\'. "
3006 << "\n#" << '\n';
3007 }
3008
3009 // define a helper to keep loops approximately dim-independent
3010 // since MathGL labels axes as x, y, z
3011 const std::string axes = "xyz";
3012
3013 // (ii) write preamble and graphing tweaks
3014 out << "\n#"
3015 << "\n# Preamble."
3016 << "\n#" << '\n';
3017
3019 out << "\nbox";
3020
3021 // deal with dimension dependent preamble; eg. default sizes and
3022 // views for MathGL (cf. gnuplot).
3023 switch (dim)
3024 {
3025 case 2:
3026 out << "\nsetsize 800 800";
3027 out << "\nrotate 0 0";
3028 break;
3029 case 3:
3030 out << "\nsetsize 800 800";
3031 out << "\nrotate 60 40";
3032 break;
3033 default:
3035 }
3036 out << '\n';
3037
3038 // (iii) write vertex ordering
3039 out << "\n#"
3040 << "\n# Vertex ordering."
3041 << "\n# list <vertex order> <vertex indices>"
3042 << "\n#" << '\n';
3043
3044 // todo: This denotes the natural ordering of vertices, but it needs
3045 // to check this is really always true for a given grid (it's not
3046 // true in @ref step_1 "step-1" grid-2 for instance).
3047 switch (dim)
3048 {
3049 case 2:
3050 out << "\nlist f 0 1 2 3" << '\n';
3051 break;
3052 case 3:
3053 out
3054 << "\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"
3055 << '\n';
3056 break;
3057 default:
3059 }
3060
3061 // (iv) write a list of vertices of cells
3062 out << "\n#"
3063 << "\n# List of vertices."
3064 << "\n# list <id> <vertices>"
3065 << "\n#" << '\n';
3066
3067 // run over all active cells and write out a list of
3068 // xyz-coordinates that correspond to vertices
3069 // No global indices in deal.II, so we make one up here.
3070 for (const auto &cell : tria.active_cell_iterators())
3071 {
3072 for (unsigned int i = 0; i < dim; ++i)
3073 {
3074 // if (cell->direction_flag ()==true)
3075 // out << "\ntrue";
3076 // else
3077 // out << "\nfalse";
3078
3079 out << "\nlist " << axes[i] << cell->active_cell_index() << " ";
3080 for (const unsigned int j : GeometryInfo<dim>::vertex_indices())
3081 out << cell->vertex(j)[i] << " ";
3082 }
3083 out << '\n';
3084 }
3085
3086 // (v) write out cells to plot as quadplot objects
3087 out << "\n#"
3088 << "\n# List of cells to quadplot."
3089 << "\n# quadplot <vertex order> <id> <style>"
3090 << "\n#" << '\n';
3091 for (unsigned int i = 0; i < tria.n_active_cells(); ++i)
3092 {
3093 out << "\nquadplot f ";
3094 for (unsigned int j = 0; j < dim; ++j)
3095 out << axes[j] << i << " ";
3096 out << "\'k#\'";
3097 }
3098 out << '\n';
3099
3100 // (vi) write footer
3101 out << "\n#"
3102 << "\n#"
3103 << "\n#" << '\n';
3104
3105 // make sure everything now gets to the output stream
3106 out.flush();
3107 AssertThrow(out.fail() == false, ExcIO());
3108}
3109
3110
3111
3112namespace
3113{
3120 template <int dim, int spacedim, typename ITERATOR, typename END>
3121 void
3122 generate_triangulation_patches(
3123 std::vector<DataOutBase::Patch<dim, spacedim>> &patches,
3124 ITERATOR cell,
3125 END end)
3126 {
3127 // convert each of the active cells into a patch
3128 for (; cell != end; ++cell)
3129 {
3131 patch.reference_cell = cell->reference_cell();
3132 patch.n_subdivisions = 1;
3133 patch.data.reinit(5, cell->n_vertices());
3134
3135 for (const unsigned int v : cell->vertex_indices())
3136 {
3137 patch.vertices[v] = cell->vertex(v);
3138 patch.data(0, v) = cell->level();
3139 patch.data(1, v) =
3140 static_cast<std::make_signed_t<types::manifold_id>>(
3141 cell->manifold_id());
3142 patch.data(2, v) =
3143 static_cast<std::make_signed_t<types::material_id>>(
3144 cell->material_id());
3145 if (cell->is_active())
3146 patch.data(3, v) =
3147 static_cast<std::make_signed_t<types::subdomain_id>>(
3148 cell->subdomain_id());
3149 else
3150 patch.data(3, v) = -1;
3151 patch.data(4, v) =
3152 static_cast<std::make_signed_t<types::subdomain_id>>(
3153 cell->level_subdomain_id());
3154 }
3155 patches.push_back(patch);
3156 }
3157 }
3158
3159
3160
3161 std::vector<std::string>
3162 triangulation_patch_data_names()
3163 {
3164 std::vector<std::string> v(5);
3165 v[0] = "level";
3166 v[1] = "manifold";
3167 v[2] = "material";
3168 v[3] = "subdomain";
3169 v[4] = "level_subdomain";
3170 return v;
3171 }
3172
3176 std::vector<typename Triangulation<3, 3>::active_line_iterator>
3177 get_boundary_edge_iterators(const Triangulation<3, 3> &tria)
3178 {
3179 std::vector<typename Triangulation<3, 3>::active_line_iterator> res;
3180
3181 std::vector<bool> flags;
3183 const_cast<Triangulation<3, 3> &>(tria).clear_user_flags_line();
3184
3185 for (auto face : tria.active_face_iterators())
3186 for (const auto l : face->line_indices())
3187 {
3188 const auto line = face->line(l);
3189 if (line->user_flag_set() || line->has_children())
3190 continue;
3191 else
3192 line->set_user_flag();
3193 if (line->at_boundary())
3194 res.emplace_back(line);
3195 }
3196 const_cast<Triangulation<3, 3> &>(tria).load_user_flags_line(flags);
3197 return res;
3198 }
3199
3200
3201
3205 template <int dim, int spacedim>
3206 std::vector<typename Triangulation<dim, spacedim>::active_line_iterator>
3207 get_boundary_edge_iterators(const Triangulation<dim, spacedim> &)
3208 {
3209 return {};
3210 }
3211
3212
3213
3218 std::vector<typename Triangulation<3, 3>::active_line_iterator>
3219 get_relevant_edge_iterators(const Triangulation<3, 3> &tria)
3220 {
3221 std::vector<typename Triangulation<3, 3>::active_line_iterator> res;
3222
3223 std::vector<bool> flags;
3225 const_cast<Triangulation<3, 3> &>(tria).clear_user_flags_line();
3226
3227 for (auto face : tria.active_face_iterators())
3228 for (const auto l : face->line_indices())
3229 {
3230 const auto line = face->line(l);
3231 if (line->user_flag_set() || line->has_children())
3232 continue;
3233 else
3234 line->set_user_flag();
3235 if (line->manifold_id() != numbers::flat_manifold_id ||
3236 (line->boundary_id() != 0 &&
3237 line->boundary_id() != numbers::invalid_boundary_id))
3238 res.emplace_back(line);
3239 }
3240 const_cast<Triangulation<3, 3> &>(tria).load_user_flags_line(flags);
3241 return res;
3242 }
3243
3244
3248 template <int dim, int spacedim>
3249 std::vector<typename Triangulation<dim, spacedim>::active_line_iterator>
3250 get_relevant_edge_iterators(const Triangulation<dim, spacedim> &)
3251 {
3252 return {};
3253 }
3254
3255
3256
3260 template <int dim, int spacedim>
3261 std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3262 get_boundary_face_iterators(const Triangulation<dim, spacedim> &tria)
3263 {
3264 std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3265 res;
3266 if (dim == 1)
3267 return res;
3268 for (auto face : tria.active_face_iterators())
3269 {
3270 if (face->boundary_id() != numbers::invalid_boundary_id)
3271 res.push_back(face);
3272 }
3273 return res;
3274 }
3275
3276
3277
3282 template <int dim, int spacedim>
3283 std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3284 get_relevant_face_iterators(const Triangulation<dim, spacedim> &tria)
3285 {
3286 std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3287 res;
3288 if (dim == 1)
3289 return res;
3290 for (auto face : tria.active_face_iterators())
3291 {
3292 if (face->manifold_id() != numbers::flat_manifold_id ||
3293 (face->boundary_id() != 0 &&
3294 face->boundary_id() != numbers::invalid_boundary_id))
3295 res.push_back(face);
3296 }
3297 return res;
3298 }
3299} // namespace
3300
3301
3302
3303template <int dim, int spacedim>
3304void
3306 std::ostream &out) const
3307{
3308 AssertThrow(out.fail() == false, ExcIO());
3309
3310 // get the positions of the vertices
3311 const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
3312
3313 const auto n_vertices = vertices.size();
3314
3315 out << "# vtk DataFile Version 3.0\n"
3316 << "Triangulation generated with deal.II\n"
3317 << "ASCII\n"
3318 << "DATASET UNSTRUCTURED_GRID\n"
3319 << "POINTS " << n_vertices << " double\n";
3320
3321 // actually write the vertices.
3322 for (const auto &v : vertices)
3323 {
3324 out << v;
3325 for (unsigned int d = spacedim + 1; d <= 3; ++d)
3326 out << " 0"; // fill with zeroes
3327 out << '\n';
3328 }
3329
3330 const auto faces = vtk_flags.output_only_relevant ?
3331 get_relevant_face_iterators(tria) :
3332 get_boundary_face_iterators(tria);
3333 const auto edges = vtk_flags.output_only_relevant ?
3334 get_relevant_edge_iterators(tria) :
3335 get_boundary_edge_iterators(tria);
3336
3338 vtk_flags.output_cells || (dim >= 2 && vtk_flags.output_faces) ||
3339 (dim >= 3 && vtk_flags.output_edges),
3340 ExcMessage(
3341 "At least one of the flags (output_cells, output_faces, output_edges) has to be enabled!"));
3342
3343 // Write cells preamble
3344 const int n_cells = (vtk_flags.output_cells ? tria.n_active_cells() : 0) +
3345 (vtk_flags.output_faces ? faces.size() : 0) +
3346 (vtk_flags.output_edges ? edges.size() : 0);
3347
3348 // VTK now expects a number telling the total storage requirement to read all
3349 // cell connectivity information. The connectivity information is read cell by
3350 // cell, first specifying how many vertices are required to describe the cell,
3351 // and then specifying the index of every vertex. This means that for every
3352 // deal.II object type, we always need n_vertices + 1 integer per cell.
3353 // Compute the total number here.
3354 int cells_size = 0;
3355
3357 for (const auto &cell : tria.active_cell_iterators())
3358 cells_size += cell->n_vertices() + 1;
3359
3361 for (const auto &face : faces)
3362 cells_size += face->n_vertices() + 1;
3363
3365 for (const auto &edge : edges)
3366 cells_size += edge->n_vertices() + 1;
3367
3368 AssertThrow(cells_size > 0, ExcMessage("No cells given to be output!"));
3369
3370 out << "\nCELLS " << n_cells << ' ' << cells_size << '\n';
3371 /*
3372 * VTK cells:
3373 *
3374 * 1 VTK_VERTEX
3375 * 3 VTK_LINE
3376 * 5 VTK_TRIANGLE
3377 * 9 VTK_QUAD
3378 * 10 VTK_TETRA
3379 * 14 VTK_PYRAMID
3380 * 13 VTK_WEDGE
3381 * 12 VTK_HEXAHEDRON
3382 *
3383 * see also: https://vtk.org/wp-content/uploads/2015/04/file-formats.pdf
3384 */
3385 static const std::array<int, 8> deal_to_vtk_cell_type = {
3386 {1, 3, 5, 9, 10, 14, 13, 12}};
3387 static const std::array<unsigned int, 8> vtk_to_deal_hypercube = {
3388 {0, 1, 3, 2, 4, 5, 7, 6}};
3389
3390 // write cells.
3392 for (const auto &cell : tria.active_cell_iterators())
3393 {
3394 out << cell->n_vertices();
3395 for (const unsigned int i : cell->vertex_indices())
3396 {
3397 out << ' ';
3398 const auto reference_cell = cell->reference_cell();
3399
3400 if ((reference_cell == ReferenceCells::Vertex) ||
3401 (reference_cell == ReferenceCells::Line) ||
3402 (reference_cell == ReferenceCells::Quadrilateral) ||
3403 (reference_cell == ReferenceCells::Hexahedron))
3404 out << cell->vertex_index(vtk_to_deal_hypercube[i]);
3405 else if ((reference_cell == ReferenceCells::Triangle) ||
3406 (reference_cell == ReferenceCells::Tetrahedron) ||
3407 (reference_cell == ReferenceCells::Wedge))
3408 out << cell->vertex_index(i);
3409 else if (reference_cell == ReferenceCells::Pyramid)
3410 {
3411 static const std::array<unsigned int, 5> permutation_table{
3412 {0, 1, 3, 2, 4}};
3413 out << cell->vertex_index(permutation_table[i]);
3414 }
3415 else
3417 }
3418 out << '\n';
3419 }
3421 for (const auto &face : faces)
3422 {
3423 out << face->n_vertices();
3424 for (const unsigned int i : face->vertex_indices())
3425 {
3426 out << ' '
3427 << face->vertex_index(GeometryInfo<dim>::vertices_per_face ==
3428 face->n_vertices() ?
3429 vtk_to_deal_hypercube[i] :
3430 i);
3431 }
3432 out << '\n';
3433 }
3435 for (const auto &edge : edges)
3436 {
3437 out << 2;
3438 for (const unsigned int i : edge->vertex_indices())
3439 out << ' ' << edge->vertex_index(i);
3440 out << '\n';
3441 }
3442
3443 // write cell types
3444 out << "\nCELL_TYPES " << n_cells << '\n';
3446 {
3447 for (const auto &cell : tria.active_cell_iterators())
3448 out << deal_to_vtk_cell_type[static_cast<int>(cell->reference_cell())]
3449 << ' ';
3450 out << '\n';
3451 }
3453 {
3454 for (const auto &face : faces)
3455 out << deal_to_vtk_cell_type[static_cast<int>(face->reference_cell())]
3456 << ' ';
3457 out << '\n';
3458 }
3460 {
3461 for (const auto &edge : edges)
3462 out << deal_to_vtk_cell_type[static_cast<int>(edge->reference_cell())]
3463 << ' ';
3464 }
3465 out << "\n\nCELL_DATA " << n_cells << '\n'
3466 << "SCALARS MaterialID int 1\n"
3467 << "LOOKUP_TABLE default\n";
3468
3469 // Now material id and boundary id
3471 {
3472 for (const auto &cell : tria.active_cell_iterators())
3473 {
3474 out << static_cast<std::make_signed_t<types::material_id>>(
3475 cell->material_id())
3476 << ' ';
3477 }
3478 out << '\n';
3479 }
3481 {
3482 for (const auto &face : faces)
3483 {
3484 out << static_cast<std::make_signed_t<types::boundary_id>>(
3485 face->boundary_id())
3486 << ' ';
3487 }
3488 out << '\n';
3489 }
3491 {
3492 for (const auto &edge : edges)
3493 {
3494 out << static_cast<std::make_signed_t<types::boundary_id>>(
3495 edge->boundary_id())
3496 << ' ';
3497 }
3498 }
3499
3500 out << "\n\nSCALARS ManifoldID int 1\n"
3501 << "LOOKUP_TABLE default\n";
3502
3503 // Now manifold id
3505 {
3506 for (const auto &cell : tria.active_cell_iterators())
3507 {
3508 out << static_cast<std::make_signed_t<types::manifold_id>>(
3509 cell->manifold_id())
3510 << ' ';
3511 }
3512 out << '\n';
3513 }
3515 {
3516 for (const auto &face : faces)
3517 {
3518 out << static_cast<std::make_signed_t<types::manifold_id>>(
3519 face->manifold_id())
3520 << ' ';
3521 }
3522 out << '\n';
3523 }
3525 {
3526 for (const auto &edge : edges)
3527 {
3528 out << static_cast<std::make_signed_t<types::manifold_id>>(
3529 edge->manifold_id())
3530 << ' ';
3531 }
3532 out << '\n';
3533 }
3534
3535 out.flush();
3536
3537 AssertThrow(out.fail() == false, ExcIO());
3538}
3539
3540
3541
3542template <int dim, int spacedim>
3543void
3545 std::ostream &out) const
3546{
3547 AssertThrow(out.fail() == false, ExcIO());
3548
3549 // convert the cells of the triangulation into a set of patches
3550 // and then have them output. since there is no data attached to
3551 // the geometry, we also do not have to provide any names, identifying
3552 // information, etc.
3553 std::vector<DataOutBase::Patch<dim, spacedim>> patches;
3554 patches.reserve(tria.n_active_cells());
3555 generate_triangulation_patches(patches, tria.begin_active(), tria.end());
3556
3559 patches,
3560 triangulation_patch_data_names(),
3561 std::vector<
3562 std::tuple<unsigned int,
3563 unsigned int,
3564 std::string,
3566 vtu_flags,
3567 out);
3569 {
3570 out << " </UnstructuredGrid>\n";
3571 out << "<dealiiData encoding=\"base64\">";
3572 std::stringstream outstring;
3573 boost::archive::binary_oarchive ia(outstring);
3574 tria.save(ia, 0);
3575 const auto compressed = Utilities::compress(outstring.str());
3576 out << Utilities::encode_base64({compressed.begin(), compressed.end()});
3577 out << "\n</dealiiData>\n";
3578 out << "</VTKFile>\n";
3579 }
3580 else
3582
3583 out << std::flush;
3584 AssertThrow(out.fail() == false, ExcIO());
3585}
3586
3587
3588
3589template <int dim, int spacedim>
3590void
3592 const Triangulation<dim, spacedim> &tria,
3593 const std::string &filename_without_extension,
3594 const bool view_levels,
3595 const bool include_artificial) const
3596{
3597 std::vector<DataOutBase::Patch<dim, spacedim>> patches;
3598 const unsigned int n_datasets = 4;
3599 std::vector<std::string> data_names;
3600 data_names.emplace_back("level");
3601 data_names.emplace_back("subdomain");
3602 data_names.emplace_back("level_subdomain");
3603 data_names.emplace_back("proc_writing");
3604
3605 const auto reference_cells = tria.get_reference_cells();
3606
3607 AssertDimension(reference_cells.size(), 1);
3608
3609 const auto &reference_cell = reference_cells[0];
3610
3611 const unsigned int n_q_points = reference_cell.n_vertices();
3612
3613 for (const auto &cell : tria.cell_iterators())
3614 {
3615 if (!view_levels)
3616 {
3617 if (cell->has_children())
3618 continue;
3619 if (!include_artificial &&
3620 cell->subdomain_id() == numbers::artificial_subdomain_id)
3621 continue;
3622 }
3623 else if (!include_artificial)
3624 {
3625 if (cell->has_children() &&
3626 cell->level_subdomain_id() == numbers::artificial_subdomain_id)
3627 continue;
3628 else if (cell->is_active() &&
3629 cell->level_subdomain_id() ==
3631 cell->subdomain_id() == numbers::artificial_subdomain_id)
3632 continue;
3633 }
3634
3636 patch.data.reinit(n_datasets, n_q_points);
3637 patch.points_are_available = false;
3638 patch.reference_cell = reference_cell;
3639
3640 for (unsigned int vertex = 0; vertex < n_q_points; ++vertex)
3641 {
3642 patch.vertices[vertex] = cell->vertex(vertex);
3643 patch.data(0, vertex) = cell->level();
3644 if (cell->is_active())
3645 patch.data(1, vertex) = static_cast<double>(
3646 static_cast<std::make_signed_t<types::subdomain_id>>(
3647 cell->subdomain_id()));
3648 else
3649 patch.data(1, vertex) = -1.0;
3650 patch.data(2, vertex) = static_cast<double>(
3651 static_cast<std::make_signed_t<types::subdomain_id>>(
3652 cell->level_subdomain_id()));
3653 patch.data(3, vertex) = tria.locally_owned_subdomain();
3654 }
3655
3656 for (auto f : reference_cell.face_indices())
3658 patches.push_back(patch);
3659 }
3660
3661 // only create .pvtu file if running in parallel
3662 // if not, just create a .vtu file with no reference
3663 // to the processor number
3664 std::string new_file = filename_without_extension + ".vtu";
3666 dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(&tria))
3667 {
3668 new_file = filename_without_extension + ".proc" +
3669 Utilities::int_to_string(tr->locally_owned_subdomain(), 4) +
3670 ".vtu";
3671
3672 // create .pvtu record
3673 if (tr->locally_owned_subdomain() == 0)
3674 {
3675 std::vector<std::string> filenames;
3676
3677 // .pvtu needs to reference the files without a relative path because
3678 // it will be written in the same directory. For this, remove any
3679 // paths from filename.
3680 std::size_t pos = filename_without_extension.find_last_of('/');
3681 if (pos == std::string::npos)
3682 pos = 0;
3683 else
3684 pos += 1;
3685 const unsigned int n_procs =
3686 Utilities::MPI::n_mpi_processes(tr->get_communicator());
3687 for (unsigned int i = 0; i < n_procs; ++i)
3688 filenames.push_back(filename_without_extension.substr(pos) +
3689 ".proc" + Utilities::int_to_string(i, 4) +
3690 ".vtu");
3691
3692 const std::string pvtu_filename =
3693 (filename_without_extension + ".pvtu");
3694 std::ofstream pvtu_output(pvtu_filename);
3695
3696 DataOut<dim, spacedim> data_out;
3697 data_out.attach_triangulation(*tr);
3698
3699 // We need a dummy vector with the names of the data values in the
3700 // .vtu files in order that the .pvtu contains reference these values
3701 Vector<float> dummy_vector(tr->n_active_cells());
3702 data_out.add_data_vector(dummy_vector, "level");
3703 data_out.add_data_vector(dummy_vector, "subdomain");
3704 data_out.add_data_vector(dummy_vector, "level_subdomain");
3705 data_out.add_data_vector(dummy_vector, "proc_writing");
3706
3707 data_out.build_patches();
3708
3709 data_out.write_pvtu_record(pvtu_output, filenames);
3710 }
3711 }
3712
3713 std::ofstream out(new_file);
3714 std::vector<
3715 std::tuple<unsigned int,
3716 unsigned int,
3717 std::string,
3719 vector_data_ranges;
3721 patches, data_names, vector_data_ranges, vtu_flags, out);
3722}
3723
3724
3725
3726unsigned int
3728{
3729 return 0;
3730}
3731
3732unsigned int
3734{
3735 return 0;
3736}
3737
3738
3739unsigned int
3741{
3742 return 0;
3743}
3744
3745unsigned int
3747{
3748 return 0;
3749}
3750
3751unsigned int
3753{
3754 return 0;
3755}
3756
3757unsigned int
3759{
3760 return 0;
3761}
3762
3763unsigned int
3765{
3766 return 0;
3767}
3768
3769unsigned int
3771{
3772 return 0;
3773}
3774
3775
3776
3777template <int dim, int spacedim>
3778unsigned int
3780{
3782 unsigned int n_faces = 0;
3783
3784 for (const auto &face : tria.active_face_iterators())
3785 if ((face->at_boundary()) && (face->boundary_id() != 0))
3786 ++n_faces;
3787
3788 return n_faces;
3789}
3790
3791
3792
3793template <int dim, int spacedim>
3794unsigned int
3796{
3797 // save the user flags for lines so
3798 // we can use these flags to track
3799 // which ones we've already counted
3800 std::vector<bool> line_flags;
3801 const_cast<::Triangulation<dim, spacedim> &>(tria).save_user_flags_line(
3802 line_flags);
3803 const_cast<::Triangulation<dim, spacedim> &>(tria)
3804 .clear_user_flags_line();
3805
3806 unsigned int n_lines = 0;
3807
3808 for (const auto &cell : tria.active_cell_iterators())
3809 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
3810 if (cell->line(l)->at_boundary() && (cell->line(l)->boundary_id() != 0) &&
3811 (cell->line(l)->user_flag_set() == false))
3812 {
3813 ++n_lines;
3814 cell->line(l)->set_user_flag();
3815 }
3816
3817 // at the end, restore the user
3818 // flags for the lines
3819 const_cast<::Triangulation<dim, spacedim> &>(tria).load_user_flags_line(
3820 line_flags);
3821
3822 return n_lines;
3823}
3824
3825
3826
3827unsigned int
3829 const unsigned int next_element_index,
3830 std::ostream &) const
3831{
3832 return next_element_index;
3833}
3834
3835
3836unsigned int
3838 const unsigned int next_element_index,
3839 std::ostream &) const
3840{
3841 return next_element_index;
3842}
3843
3844unsigned int
3846 const unsigned int next_element_index,
3847 std::ostream &) const
3848{
3849 return next_element_index;
3850}
3851
3852
3853unsigned int
3855 const unsigned int next_element_index,
3856 std::ostream &) const
3857{
3858 return next_element_index;
3859}
3860
3861unsigned int
3863 const unsigned int next_element_index,
3864 std::ostream &) const
3865{
3866 return next_element_index;
3867}
3868
3869
3870unsigned int
3872 const unsigned int next_element_index,
3873 std::ostream &) const
3874{
3875 return next_element_index;
3876}
3877
3878
3879unsigned int
3881 const unsigned int next_element_index,
3882 std::ostream &) const
3883{
3884 return next_element_index;
3885}
3886
3887unsigned int
3889 const unsigned int next_element_index,
3890 std::ostream &) const
3891{
3892 return next_element_index;
3893}
3894
3895
3896
3897template <int dim, int spacedim>
3898unsigned int
3900 const unsigned int next_element_index,
3901 std::ostream &out) const
3902{
3903 unsigned int current_element_index = next_element_index;
3904
3905 for (const auto &face : tria.active_face_iterators())
3906 if (face->at_boundary() && (face->boundary_id() != 0))
3907 {
3908 out << current_element_index << ' '
3909 << face->reference_cell().gmsh_element_type() << ' ';
3910 out << static_cast<unsigned int>(face->boundary_id()) << ' '
3911 << static_cast<unsigned int>(face->boundary_id()) << ' '
3912 << face->n_vertices();
3913 // note: vertex numbers are 1-base
3914 for (const unsigned int vertex : face->vertex_indices())
3915 {
3916 if (face->reference_cell() == ReferenceCells::Quadrilateral)
3917 out << ' '
3918 << face->vertex_index(
3920 1;
3921 else if ((face->reference_cell() == ReferenceCells::Triangle) ||
3922 (face->reference_cell() == ReferenceCells::Line))
3923 out << ' ' << face->vertex_index(vertex) + 1;
3924 else
3926 }
3927 out << '\n';
3928
3929 ++current_element_index;
3930 }
3931 return current_element_index;
3932}
3933
3934
3935
3936template <int dim, int spacedim>
3937unsigned int
3939 const unsigned int next_element_index,
3940 std::ostream &out) const
3941{
3942 unsigned int current_element_index = next_element_index;
3943 // save the user flags for lines so
3944 // we can use these flags to track
3945 // which ones we've already taken
3946 // care of
3947 std::vector<bool> line_flags;
3948 const_cast<::Triangulation<dim, spacedim> &>(tria).save_user_flags_line(
3949 line_flags);
3950 const_cast<::Triangulation<dim, spacedim> &>(tria)
3951 .clear_user_flags_line();
3952
3953 for (const auto &cell : tria.active_cell_iterators())
3954 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
3955 if (cell->line(l)->at_boundary() && (cell->line(l)->boundary_id() != 0) &&
3956 (cell->line(l)->user_flag_set() == false))
3957 {
3958 out << next_element_index << ' '
3960 out << static_cast<unsigned int>(cell->line(l)->boundary_id()) << ' '
3961 << static_cast<unsigned int>(cell->line(l)->boundary_id())
3962 << " 2 "; // two vertex indices to follow
3963 // note: vertex numbers are 1-base
3964 for (unsigned int vertex = 0; vertex < 2; ++vertex)
3965 out << ' '
3966 << cell->line(l)->vertex_index(
3968 1;
3969 out << '\n';
3970
3971 // move on to the next line
3972 // but mark the current one
3973 // as taken care of
3974 ++current_element_index;
3975 cell->line(l)->set_user_flag();
3976 }
3977
3978 // at the end, restore the user
3979 // flags for the lines
3980 const_cast<::Triangulation<dim, spacedim> &>(tria).load_user_flags_line(
3981 line_flags);
3982
3983 return current_element_index;
3984}
3985
3986
3987
3988unsigned int
3990 const unsigned int next_element_index,
3991 std::ostream &) const
3992{
3993 return next_element_index;
3994}
3995
3996unsigned int
3998 const unsigned int next_element_index,
3999 std::ostream &) const
4000{
4001 return next_element_index;
4002}
4003
4004unsigned int
4006 const unsigned int next_element_index,
4007 std::ostream &) const
4008{
4009 return next_element_index;
4010}
4011
4012unsigned int
4014 const unsigned int next_element_index,
4015 std::ostream &) const
4016{
4017 return next_element_index;
4018}
4019
4020unsigned int
4022 const unsigned int next_element_index,
4023 std::ostream &) const
4024{
4025 return next_element_index;
4026}
4027
4028
4029unsigned int
4031 const unsigned int next_element_index,
4032 std::ostream &) const
4033{
4034 return next_element_index;
4035}
4036
4037
4038unsigned int
4040 const unsigned int next_element_index,
4041 std::ostream &) const
4042{
4043 return next_element_index;
4044}
4045
4046unsigned int
4048 const unsigned int next_element_index,
4049 std::ostream &) const
4050{
4051 return next_element_index;
4052}
4053
4054
4055
4056template <int dim, int spacedim>
4057unsigned int
4059 const unsigned int next_element_index,
4060 std::ostream &out) const
4061{
4062 unsigned int current_element_index = next_element_index;
4064
4065 for (const auto &face : tria.active_face_iterators())
4066 if (face->at_boundary() && (face->boundary_id() != 0))
4067 {
4068 out << current_element_index << " "
4069 << static_cast<unsigned int>(face->boundary_id()) << " ";
4070 switch (dim)
4071 {
4072 case 2:
4073 out << "line ";
4074 break;
4075 case 3:
4076 out << "quad ";
4077 break;
4078 default:
4080 }
4081 // note: vertex numbers are 1-base
4082 for (unsigned int vertex = 0;
4083 vertex < GeometryInfo<dim>::vertices_per_face;
4084 ++vertex)
4085 out << face->vertex_index(
4087 1
4088 << ' ';
4089 out << '\n';
4090
4091 ++current_element_index;
4092 }
4093 return current_element_index;
4094}
4095
4096
4097
4098template <int dim, int spacedim>
4099unsigned int
4101 const unsigned int next_element_index,
4102 std::ostream &out) const
4103{
4104 unsigned int current_element_index = next_element_index;
4105 // save the user flags for lines so
4106 // we can use these flags to track
4107 // which ones we've already taken
4108 // care of
4109 std::vector<bool> line_flags;
4110 const_cast<::Triangulation<dim, spacedim> &>(tria).save_user_flags_line(
4111 line_flags);
4112 const_cast<::Triangulation<dim, spacedim> &>(tria)
4113 .clear_user_flags_line();
4114
4115 for (const auto &cell : tria.active_cell_iterators())
4116 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
4117 if (cell->line(l)->at_boundary() && (cell->line(l)->boundary_id() != 0) &&
4118 (cell->line(l)->user_flag_set() == false))
4119 {
4120 out << current_element_index << " "
4121 << static_cast<unsigned int>(cell->line(l)->boundary_id())
4122 << " line ";
4123 // note: vertex numbers in ucd format are 1-base
4124 for (unsigned int vertex = 0; vertex < 2; ++vertex)
4125 out << cell->line(l)->vertex_index(
4127 1
4128 << ' ';
4129 out << '\n';
4130
4131 // move on to the next line
4132 // but mark the current one
4133 // as taken care of
4134 ++current_element_index;
4135 cell->line(l)->set_user_flag();
4136 }
4137
4138 // at the end, restore the user
4139 // flags for the lines
4140 const_cast<::Triangulation<dim, spacedim> &>(tria).load_user_flags_line(
4141 line_flags);
4142 return current_element_index;
4143}
4144
4145
4146
4147namespace internal
4148{
4149 namespace
4150 {
4159 template <int spacedim>
4160 void
4161 remove_colinear_points(std::vector<Point<spacedim>> &points)
4162 {
4163 while (points.size() > 2)
4164 {
4165 Tensor<1, spacedim> first_difference = points[1] - points[0];
4166 first_difference /= first_difference.norm();
4167 Tensor<1, spacedim> second_difference = points[2] - points[1];
4168 second_difference /= second_difference.norm();
4169 // If the three points are colinear then remove the middle one.
4170 if ((first_difference - second_difference).norm() < 1e-10)
4171 points.erase(points.begin() + 1);
4172 else
4173 break;
4174 }
4175 }
4176
4177
4178
4179 template <int spacedim>
4180 void
4181 write_gnuplot(const ::Triangulation<1, spacedim> &tria,
4182 std::ostream &out,
4183 const Mapping<1, spacedim> *,
4184 const GridOutFlags::Gnuplot &gnuplot_flags)
4185 {
4186 AssertThrow(out.fail() == false, ExcIO());
4187
4188 for (const auto &cell : tria.active_cell_iterators())
4189 {
4190 if (gnuplot_flags.write_cell_numbers)
4191 out << "# cell " << cell << '\n';
4192
4193 out << cell->vertex(0) << ' ' << cell->level() << ' '
4194 << cell->material_id() << '\n'
4195 << cell->vertex(1) << ' ' << cell->level() << ' '
4196 << cell->material_id() << '\n'
4197 << "\n\n";
4198 }
4199
4200 // make sure everything now gets to
4201 // disk
4202 out.flush();
4203
4204 AssertThrow(out.fail() == false, ExcIO());
4205 }
4206
4207
4208
4209 template <int spacedim>
4210 void
4211 write_gnuplot(const ::Triangulation<2, spacedim> &tria,
4212 std::ostream &out,
4213 const Mapping<2, spacedim> *mapping,
4214 const GridOutFlags::Gnuplot &gnuplot_flags)
4215 {
4216 AssertThrow(out.fail() == false, ExcIO());
4217
4218 const int dim = 2;
4219
4220 const unsigned int n_additional_points =
4221 gnuplot_flags.n_extra_curved_line_points;
4222 const unsigned int n_points = 2 + n_additional_points;
4223
4224 // If we need to plot curved lines then generate a quadrature formula to
4225 // place points via the mapping
4226 Quadrature<dim> q_projector;
4227 std::vector<Point<dim - 1>> boundary_points;
4228 if (mapping != nullptr)
4229 {
4230 boundary_points.resize(n_points);
4231 boundary_points[0][0] = 0;
4232 boundary_points[n_points - 1][0] = 1;
4233 for (unsigned int i = 1; i < n_points - 1; ++i)
4234 boundary_points[i][0] = 1. * i / (n_points - 1);
4235
4236 std::vector<double> dummy_weights(n_points, 1. / n_points);
4237 Quadrature<dim - 1> quadrature(boundary_points, dummy_weights);
4238
4240 ReferenceCells::get_hypercube<dim>(), quadrature);
4241 }
4242
4243 static constexpr std::array<unsigned int, 8> local_vertex_numbering = {
4244 {0, 1, 5, 4, 2, 3, 7, 6}};
4245 for (const auto &cell : tria.active_cell_iterators())
4246 {
4247 if (gnuplot_flags.write_cell_numbers)
4248 out << "# cell " << cell << '\n';
4249
4250 if (mapping == nullptr ||
4251 (dim == spacedim ?
4252 (!cell->at_boundary() && !gnuplot_flags.curved_inner_cells) :
4253 // ignore checking for boundary or interior cells in the codim
4254 // 1 case: 'or false' is a no-op
4255 false))
4256 {
4257 // write out the four sides of this cell by putting the four
4258 // points (+ the initial point again) in a row and lifting the
4259 // drawing pencil at the end
4260 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
4261 out << cell->vertex(dim == 3 ?
4262 local_vertex_numbering[i] :
4263 GeometryInfo<dim>::ucd_to_deal[i])
4264 << ' ' << cell->level() << ' ' << cell->material_id()
4265 << '\n';
4266 out << cell->vertex(0) << ' ' << cell->level() << ' '
4267 << cell->material_id() << '\n'
4268 << '\n' // double new line for gnuplot 3d plots
4269 << '\n';
4270 }
4271 else
4272 // cell is at boundary and we are to treat curved boundaries. so
4273 // loop over all faces and draw them as small pieces of lines
4274 {
4275 for (const unsigned int face_no :
4276 GeometryInfo<dim>::face_indices())
4277 {
4278 const typename ::Triangulation<dim,
4279 spacedim>::face_iterator
4280 face = cell->face(face_no);
4281 if (dim != spacedim || face->at_boundary() ||
4282 gnuplot_flags.curved_inner_cells)
4283 {
4284 // Save the points on each face to a vector and then try
4285 // to remove colinear points that won't show up in the
4286 // generated plot.
4287 std::vector<Point<spacedim>> line_points;
4288 // compute offset of quadrature points within set of
4289 // projected points
4290 const unsigned int offset = face_no * n_points;
4291 for (unsigned int i = 0; i < n_points; ++i)
4292 line_points.push_back(
4293 mapping->transform_unit_to_real_cell(
4294 cell, q_projector.point(offset + i)));
4295 internal::remove_colinear_points(line_points);
4296
4297 for (const Point<spacedim> &point : line_points)
4298 out << point << ' ' << cell->level() << ' '
4299 << cell->material_id() << '\n';
4300
4301 out << '\n' << '\n';
4302 }
4303 else
4304 {
4305 // if, however, the face is not at the boundary and we
4306 // don't want to curve anything, then draw it as usual
4307 out << face->vertex(0) << ' ' << cell->level() << ' '
4308 << cell->material_id() << '\n'
4309 << face->vertex(1) << ' ' << cell->level() << ' '
4310 << cell->material_id() << '\n'
4311 << '\n'
4312 << '\n';
4313 }
4314 }
4315 }
4316 }
4317
4318 // make sure everything now gets to disk
4319 out.flush();
4320
4321 AssertThrow(out.fail() == false, ExcIO());
4322 }
4323
4324
4325
4326 template <int spacedim>
4327 void
4328 write_gnuplot(const ::Triangulation<3, spacedim> &tria,
4329 std::ostream &out,
4330 const Mapping<3, spacedim> *mapping,
4331 const GridOutFlags::Gnuplot &gnuplot_flags)
4332 {
4333 AssertThrow(out.fail() == false, ExcIO());
4334
4335 const int dim = 3;
4336
4337 const unsigned int n_additional_points =
4338 gnuplot_flags.n_extra_curved_line_points;
4339 const unsigned int n_points = 2 + n_additional_points;
4340
4341 // If we need to plot curved lines then generate a quadrature formula to
4342 // place points via the mapping
4343 std::unique_ptr<Quadrature<dim>> q_projector;
4344 std::vector<Point<1>> boundary_points;
4345 if (mapping != nullptr)
4346 {
4347 boundary_points.resize(n_points);
4348 boundary_points[0][0] = 0;
4349 boundary_points[n_points - 1][0] = 1;
4350 for (unsigned int i = 1; i < n_points - 1; ++i)
4351 boundary_points[i][0] = 1. * i / (n_points - 1);
4352
4353 std::vector<double> dummy_weights(n_points, 1. / n_points);
4354 Quadrature<1> quadrature1d(boundary_points, dummy_weights);
4355
4356 // tensor product of points, only one copy
4357 QIterated<dim - 1> quadrature(quadrature1d, 1);
4358 q_projector = std::make_unique<Quadrature<dim>>(
4360 ReferenceCells::get_hypercube<dim>(), quadrature));
4361 }
4362
4363 for (const auto &cell : tria.active_cell_iterators())
4364 {
4365 if (gnuplot_flags.write_cell_numbers)
4366 out << "# cell " << cell << '\n';
4367
4368 if (mapping == nullptr || n_points == 2 ||
4369 (!cell->has_boundary_lines() &&
4370 !gnuplot_flags.curved_inner_cells))
4371 {
4372 if (cell->reference_cell() == ReferenceCells::Hexahedron)
4373 {
4374 // front face
4375 out << cell->vertex(0) << ' ' << cell->level() << ' '
4376 << cell->material_id() << '\n'
4377 << cell->vertex(1) << ' ' << cell->level() << ' '
4378 << cell->material_id() << '\n'
4379 << cell->vertex(5) << ' ' << cell->level() << ' '
4380 << cell->material_id() << '\n'
4381 << cell->vertex(4) << ' ' << cell->level() << ' '
4382 << cell->material_id() << '\n'
4383 << cell->vertex(0) << ' ' << cell->level() << ' '
4384 << cell->material_id() << '\n'
4385 << '\n';
4386 // back face
4387 out << cell->vertex(2) << ' ' << cell->level() << ' '
4388 << cell->material_id() << '\n'
4389 << cell->vertex(3) << ' ' << cell->level() << ' '
4390 << cell->material_id() << '\n'
4391 << cell->vertex(7) << ' ' << cell->level() << ' '
4392 << cell->material_id() << '\n'
4393 << cell->vertex(6) << ' ' << cell->level() << ' '
4394 << cell->material_id() << '\n'
4395 << cell->vertex(2) << ' ' << cell->level() << ' '
4396 << cell->material_id() << '\n'
4397 << '\n';
4398
4399 // now for the four connecting lines
4400 out << cell->vertex(0) << ' ' << cell->level() << ' '
4401 << cell->material_id() << '\n'
4402 << cell->vertex(2) << ' ' << cell->level() << ' '
4403 << cell->material_id() << '\n'
4404 << '\n';
4405 out << cell->vertex(1) << ' ' << cell->level() << ' '
4406 << cell->material_id() << '\n'
4407 << cell->vertex(3) << ' ' << cell->level() << ' '
4408 << cell->material_id() << '\n'
4409 << '\n';
4410 out << cell->vertex(5) << ' ' << cell->level() << ' '
4411 << cell->material_id() << '\n'
4412 << cell->vertex(7) << ' ' << cell->level() << ' '
4413 << cell->material_id() << '\n'
4414 << '\n';
4415 out << cell->vertex(4) << ' ' << cell->level() << ' '
4416 << cell->material_id() << '\n'
4417 << cell->vertex(6) << ' ' << cell->level() << ' '
4418 << cell->material_id() << '\n'
4419 << '\n';
4420 }
4421 else if (cell->reference_cell() == ReferenceCells::Tetrahedron)
4422 {
4423 // Draw the tetrahedron as a two collections of lines.
4424 for (const unsigned int v : {0, 1, 2, 0, 3, 2})
4425 {
4426 out << cell->vertex(v) << ' ' << cell->level() << ' '
4427 << cell->material_id() << '\n';
4428 }
4429 out << '\n'; // end of first line
4430
4431 for (const unsigned int v : {3, 1})
4432 {
4433 out << cell->vertex(v) << ' ' << cell->level() << ' '
4434 << cell->material_id() << '\n';
4435 }
4436 out << '\n'; // end of second line
4437 }
4438 else if (cell->reference_cell() == ReferenceCells::Wedge)
4439 {
4440 // Draw the wedge as a collection of three
4441 // lines. The first one wraps around the base,
4442 // goes up to the top, and wraps around that. The
4443 // second and third are just individual lines
4444 // going from base to top.
4445 for (const unsigned int v : {0, 1, 2, 0, 3, 4, 5, 3})
4446 {
4447 out << cell->vertex(v) << ' ' << cell->level() << ' '
4448 << cell->material_id() << '\n';
4449 }
4450 out << '\n'; // end of first line
4451
4452 for (const unsigned int v : {1, 4})
4453 {
4454 out << cell->vertex(v) << ' ' << cell->level() << ' '
4455 << cell->material_id() << '\n';
4456 }
4457 out << '\n'; // end of second line
4458
4459 for (const unsigned int v : {2, 5})
4460 {
4461 out << cell->vertex(v) << ' ' << cell->level() << ' '
4462 << cell->material_id() << '\n';
4463 }
4464 out << '\n'; // end of third line
4465 }
4466 else if (cell->reference_cell() == ReferenceCells::Pyramid)
4467 {
4468 // Draw the pyramid as a collections of two lines.
4469 for (const unsigned int v : {0, 1, 3, 2, 0, 4, 1})
4470 {
4471 out << cell->vertex(v) << ' ' << cell->level() << ' '
4472 << cell->material_id() << '\n';
4473 }
4474 out << '\n'; // end of first line
4475
4476 for (const unsigned int v : {2, 4, 3})
4477 {
4478 out << cell->vertex(v) << ' ' << cell->level() << ' '
4479 << cell->material_id() << '\n';
4480 }
4481 out << '\n'; // end of second line
4482 }
4483 else
4485 }
4486 else // need to handle curved boundaries
4487 {
4488 Assert(cell->reference_cell() == ReferenceCells::Hexahedron,
4490 for (const unsigned int face_no :
4491 GeometryInfo<dim>::face_indices())
4492 {
4493 const typename ::Triangulation<dim,
4494 spacedim>::face_iterator
4495 face = cell->face(face_no);
4496
4497 if (face->at_boundary() &&
4498 gnuplot_flags.write_additional_boundary_lines)
4499 {
4500 const unsigned int offset = face_no * n_points * n_points;
4501 for (unsigned int i = 0; i < n_points - 1; ++i)
4502 for (unsigned int j = 0; j < n_points - 1; ++j)
4503 {
4504 const Point<spacedim> p0 =
4505 mapping->transform_unit_to_real_cell(
4506 cell,
4507 q_projector->point(offset + i * n_points + j));
4508 out << p0 << ' ' << cell->level() << ' '
4509 << cell->material_id() << '\n';
4510 out << (mapping->transform_unit_to_real_cell(
4511 cell,
4512 q_projector->point(
4513 offset + (i + 1) * n_points + j)))
4514 << ' ' << cell->level() << ' '
4515 << cell->material_id() << '\n';
4516 out << (mapping->transform_unit_to_real_cell(
4517 cell,
4518 q_projector->point(
4519 offset + (i + 1) * n_points + j + 1)))
4520 << ' ' << cell->level() << ' '
4521 << cell->material_id() << '\n';
4522 out << (mapping->transform_unit_to_real_cell(
4523 cell,
4524 q_projector->point(offset + i * n_points +
4525 j + 1)))
4526 << ' ' << cell->level() << ' '
4527 << cell->material_id() << '\n';
4528 // and the first point again
4529 out << p0 << ' ' << cell->level() << ' '
4530 << cell->material_id() << '\n';
4531 out << '\n' << '\n';
4532 }
4533 }
4534 else
4535 {
4536 for (unsigned int l = 0;
4537 l < GeometryInfo<dim>::lines_per_face;
4538 ++l)
4539 {
4540 const typename ::Triangulation<dim, spacedim>::
4541 line_iterator line = face->line(l);
4542
4543 const Point<spacedim> &v0 = line->vertex(0),
4544 &v1 = line->vertex(1);
4545 if (line->at_boundary() ||
4546 gnuplot_flags.curved_inner_cells)
4547 {
4548 // Save the points on each face to a vector and
4549 // then try to remove colinear points that won't
4550 // show up in the generated plot.
4551 std::vector<Point<spacedim>> line_points;
4552 // transform_real_to_unit_cell could be replaced
4553 // by using QProjector<dim>::project_to_line
4554 // which is not yet implemented
4555 const Point<spacedim>
4556 u0 = mapping->transform_real_to_unit_cell(cell,
4557 v0),
4558 u1 = mapping->transform_real_to_unit_cell(cell,
4559 v1);
4560 for (unsigned int i = 0; i < n_points; ++i)
4561 line_points.push_back(
4562 mapping->transform_unit_to_real_cell(
4563 cell,
4564 (1 - boundary_points[i][0]) * u0 +
4565 boundary_points[i][0] * u1));
4566 internal::remove_colinear_points(line_points);
4567 for (const Point<spacedim> &point : line_points)
4568 out << point << ' ' << cell->level() << ' '
4569 << static_cast<unsigned int>(
4570 cell->material_id())
4571 << '\n';
4572 }
4573 else
4574 out << v0 << ' ' << cell->level() << ' '
4575 << cell->material_id() << '\n'
4576 << v1 << ' ' << cell->level() << ' '
4577 << cell->material_id() << '\n';
4578
4579 out << '\n' << '\n';
4580 }
4581 }
4582 }
4583 }
4584 }
4585
4586 // make sure everything now gets to disk
4587 out.flush();
4588
4589 AssertThrow(out.fail() == false, ExcIO());
4590 }
4591 } // namespace
4592} // namespace internal
4593
4594
4595
4596template <int dim, int spacedim>
4597void
4599 std::ostream &out,
4600 const Mapping<dim, spacedim> *mapping) const
4601{
4602 internal::write_gnuplot(tria, out, mapping, gnuplot_flags);
4603}
4604
4605
4606
4607namespace internal
4608{
4609 namespace
4610 {
4611 struct LineEntry
4612 {
4616 unsigned int level;
4617 LineEntry(const Point<2> &f,
4618 const Point<2> &s,
4619 const bool c,
4620 const unsigned int l)
4621 : first(f)
4622 , second(s)
4623 , colorize(c)
4624 , level(l)
4625 {}
4626 };
4627
4628
4629 void
4630 write_eps(const ::Triangulation<1> &,
4631 std::ostream &,
4632 const Mapping<1> *,
4633 const GridOutFlags::Eps<2> &,
4634 const GridOutFlags::Eps<3> &)
4635 {
4637 }
4638
4639 void
4640 write_eps(const ::Triangulation<1, 2> &,
4641 std::ostream &,
4642 const Mapping<1, 2> *,
4643 const GridOutFlags::Eps<2> &,
4644 const GridOutFlags::Eps<3> &)
4645 {
4647 }
4648
4649 void
4650 write_eps(const ::Triangulation<1, 3> &,
4651 std::ostream &,
4652 const Mapping<1, 3> *,
4653 const GridOutFlags::Eps<2> &,
4654 const GridOutFlags::Eps<3> &)
4655 {
4657 }
4658
4659 void
4660 write_eps(const ::Triangulation<2, 3> &,
4661 std::ostream &,
4662 const Mapping<2, 3> *,
4663 const GridOutFlags::Eps<2> &,
4664 const GridOutFlags::Eps<3> &)
4665 {
4667 }
4668
4669
4670
4671 template <int dim, int spacedim>
4672 void
4673 write_eps(const ::Triangulation<dim, spacedim> &tria,
4674 std::ostream &out,
4675 const Mapping<dim, spacedim> *mapping,
4676 const GridOutFlags::Eps<2> &eps_flags_2,
4677 const GridOutFlags::Eps<3> &eps_flags_3)
4678 {
4679 using LineList = std::list<LineEntry>;
4680
4681 // We should never get here in 1d since this function is overloaded for
4682 // all dim == 1 cases.
4683 Assert(dim == 2 || dim == 3, ExcInternalError());
4684
4685 // Copy, with an object slice, something containing the flags common to
4686 // all dimensions in order to avoid the recurring distinctions between
4687 // the different eps_flags present.
4688 const GridOutFlags::EpsFlagsBase eps_flags_base =
4689 dim == 2 ?
4690 static_cast<const GridOutFlags::EpsFlagsBase &>(eps_flags_2) :
4691 static_cast<const GridOutFlags::EpsFlagsBase &>(eps_flags_3);
4692
4693 AssertThrow(out.fail() == false, ExcIO());
4694 const unsigned int n_points = eps_flags_base.n_boundary_face_points;
4695
4696 // make up a list of lines by which
4697 // we will construct the triangulation
4698 //
4699 // this part unfortunately is a bit
4700 // dimension dependent, so we have to
4701 // treat every dimension different.
4702 // however, by directly producing
4703 // the lines to be printed, i.e. their
4704 // 2d images, we can later do the
4705 // actual output dimension independent
4706 // again
4707 LineList line_list;
4708
4709 switch (dim)
4710 {
4711 case 1:
4712 {
4714 break;
4715 }
4716
4717 case 2:
4718 {
4719 for (const auto &cell : tria.active_cell_iterators())
4720 for (const unsigned int line_no : cell->line_indices())
4721 {
4722 typename ::Triangulation<dim, spacedim>::line_iterator
4723 line = cell->line(line_no);
4724
4725 // first treat all
4726 // interior lines and
4727 // make up a list of
4728 // them. if curved
4729 // lines shall not be
4730 // supported (i.e. no
4731 // mapping is
4732 // provided), then also
4733 // treat all other
4734 // lines
4735 if (!line->has_children() &&
4736 (mapping == nullptr || !line->at_boundary()))
4737 // one would expect
4738 // make_pair(line->vertex(0),
4739 // line->vertex(1))
4740 // here, but that is
4741 // not dimension
4742 // independent, since
4743 // vertex(i) is
4744 // Point<dim>, but we
4745 // want a Point<2>.
4746 // in fact, whenever
4747 // we're here, the
4748 // vertex is a
4749 // Point<dim>, but
4750 // the compiler does
4751 // not know
4752 // this. hopefully,
4753 // the compiler will
4754 // optimize away this
4755 // little kludge
4756 line_list.emplace_back(
4757 Point<2>(line->vertex(0)[0], line->vertex(0)[1]),
4758 Point<2>(line->vertex(1)[0], line->vertex(1)[1]),
4759 line->user_flag_set(),
4760 cell->level());
4761 }
4762
4763 // next if we are to treat
4764 // curved boundaries
4765 // specially, then add lines
4766 // to the list consisting of
4767 // pieces of the boundary
4768 // lines
4769 if (mapping != nullptr)
4770 {
4771 // to do so, first
4772 // generate a sequence of
4773 // points on a face and
4774 // project them onto the
4775 // faces of a unit cell
4776 std::vector<Point<dim - 1>> boundary_points(n_points);
4777
4778 for (unsigned int i = 0; i < n_points; ++i)
4779 boundary_points[i][0] = 1. * (i + 1) / (n_points + 1);
4780
4781 Quadrature<dim - 1> quadrature(boundary_points);
4782 Quadrature<dim> q_projector(
4784 ReferenceCells::get_hypercube<dim>(), quadrature));
4785
4786 // next loop over all
4787 // boundary faces and
4788 // generate the info from
4789 // them
4790 for (const auto &cell : tria.active_cell_iterators())
4791 for (const unsigned int face_no :
4792 GeometryInfo<dim>::face_indices())
4793 {
4794 const typename ::Triangulation<dim, spacedim>::
4795 face_iterator face = cell->face(face_no);
4796
4797 if (face->at_boundary())
4798 {
4799 Point<dim> p0_dim(face->vertex(0));
4800 Point<2> p0(p0_dim[0], p0_dim[1]);
4801
4802 // loop over
4803 // all pieces
4804 // of the line
4805 // and generate
4806 // line-lets
4807 const unsigned int offset = face_no * n_points;
4808 for (unsigned int i = 0; i < n_points; ++i)
4809 {
4810 const Point<dim> p1_dim(
4811 mapping->transform_unit_to_real_cell(
4812 cell, q_projector.point(offset + i)));
4813 const Point<2> p1(p1_dim[0], p1_dim[1]);
4814
4815 line_list.emplace_back(p0,
4816 p1,
4817 face->user_flag_set(),
4818 cell->level());
4819 p0 = p1;
4820 }
4821
4822 // generate last piece
4823 const Point<dim> p1_dim(face->vertex(1));
4824 const Point<2> p1(p1_dim[0], p1_dim[1]);
4825 line_list.emplace_back(p0,
4826 p1,
4827 face->user_flag_set(),
4828 cell->level());
4829 }
4830 }
4831 }
4832
4833 break;
4834 }
4835
4836 case 3:
4837 {
4838 // curved boundary output
4839 // presently not supported
4840 Assert(mapping == nullptr, ExcNotImplemented());
4841
4842 // loop over all lines and compute their
4843 // projection on the plane perpendicular
4844 // to the direction of sight
4845
4846 // direction of view equals the unit
4847 // vector of the position of the
4848 // spectator to the origin.
4849 //
4850 // we chose here the viewpoint as in
4851 // gnuplot as default.
4852 //
4853 // TODO:[WB] Fix a potential problem with viewing angles in 3d Eps
4854 // GridOut
4855 // note: the following might be wrong
4856 // if one of the base vectors below
4857 // is in direction of the viewer, but
4858 // I am too tired at present to fix
4859 // this
4860 const double pi = numbers::PI;
4861 const double z_angle = eps_flags_3.azimut_angle;
4862 const double turn_angle = eps_flags_3.turn_angle;
4863 const Point<dim> view_direction(
4864 -std::sin(z_angle * 2. * pi / 360.) *
4865 std::sin(turn_angle * 2. * pi / 360.),
4866 +std::sin(z_angle * 2. * pi / 360.) *
4867 std::cos(turn_angle * 2. * pi / 360.),
4868 -std::cos(z_angle * 2. * pi / 360.));
4869
4870 // decide about the two unit vectors
4871 // in this plane. we chose the first one
4872 // to be the projection of the z-axis
4873 // to this plane
4874 const Tensor<1, dim> vector1 =
4875 Point<dim>(0, 0, 1) -
4876 ((Point<dim>(0, 0, 1) * view_direction) * view_direction);
4877 const Tensor<1, dim> unit_vector1 = vector1 / vector1.norm();
4878
4879 // now the third vector is fixed. we
4880 // chose the projection of a more or
4881 // less arbitrary vector to the plane
4882 // perpendicular to the first one
4883 const Tensor<1, dim> vector2 =
4884 (Point<dim>(1, 0, 0) -
4885 ((Point<dim>(1, 0, 0) * view_direction) * view_direction) -
4886 ((Point<dim>(1, 0, 0) * unit_vector1) * unit_vector1));
4887 const Tensor<1, dim> unit_vector2 = vector2 / vector2.norm();
4888
4889
4890 for (const auto &cell : tria.active_cell_iterators())
4891 for (const unsigned int line_no : cell->line_indices())
4892 {
4893 typename ::Triangulation<dim, spacedim>::line_iterator
4894 line = cell->line(line_no);
4895 line_list.emplace_back(
4896 Point<2>(line->vertex(0) * unit_vector2,
4897 line->vertex(0) * unit_vector1),
4898 Point<2>(line->vertex(1) * unit_vector2,
4899 line->vertex(1) * unit_vector1),
4900 line->user_flag_set(),
4901 cell->level());
4902 }
4903
4904 break;
4905 }
4906
4907 default:
4909 }
4910
4911
4912
4913 // find out minimum and maximum x and
4914 // y coordinates to compute offsets
4915 // and scaling factors
4916 double x_min = tria.begin_active()->vertex(0)[0];
4917 double x_max = x_min;
4918 double y_min = tria.begin_active()->vertex(0)[1];
4919 double y_max = y_min;
4920 unsigned int max_level = line_list.begin()->level;
4921
4922 for (LineList::const_iterator line = line_list.begin();
4923 line != line_list.end();
4924 ++line)
4925 {
4926 x_min = std::min(x_min, line->first[0]);
4927 x_min = std::min(x_min, line->second[0]);
4928
4929 x_max = std::max(x_max, line->first[0]);
4930 x_max = std::max(x_max, line->second[0]);
4931
4932 y_min = std::min(y_min, line->first[1]);
4933 y_min = std::min(y_min, line->second[1]);
4934
4935 y_max = std::max(y_max, line->first[1]);
4936 y_max = std::max(y_max, line->second[1]);
4937
4938 max_level = std::max(max_level, line->level);
4939 }
4940
4941 // scale in x-direction such that
4942 // in the output 0 <= x <= 300.
4943 // don't scale in y-direction to
4944 // preserve the shape of the
4945 // triangulation
4946 const double scale =
4947 (eps_flags_base.size /
4948 (eps_flags_base.size_type == GridOutFlags::EpsFlagsBase::width ?
4949 x_max - x_min :
4950 y_min - y_max));
4951
4952
4953 // now write preamble
4954 {
4955 // block this to have local
4956 // variables destroyed after
4957 // use
4958 std::time_t time1 = std::time(nullptr);
4959 std::tm *time = std::localtime(&time1);
4960 out << "%!PS-Adobe-2.0 EPSF-1.2" << '\n'
4961 << "%%Title: deal.II Output" << '\n'
4962 << "%%Creator: the deal.II library" << '\n'
4963 << "%%Creation Date: " << time->tm_year + 1900 << "/"
4964 << time->tm_mon + 1 << "/" << time->tm_mday << " - "
4965 << time->tm_hour << ":" << std::setw(2) << time->tm_min << ":"
4966 << std::setw(2) << time->tm_sec << '\n'
4967 << "%%BoundingBox: "
4968 // lower left corner
4969 << "0 0 "
4970 // upper right corner
4971 << static_cast<unsigned int>(
4972 std::floor(((x_max - x_min) * scale) + 1))
4973 << ' '
4974 << static_cast<unsigned int>(
4975 std::floor(((y_max - y_min) * scale) + 1))
4976 << '\n';
4977
4978 // define some abbreviations to keep
4979 // the output small:
4980 // m=move turtle to
4981 // x=execute line stroke
4982 // b=black pen
4983 // r=red pen
4984 out << "/m {moveto} bind def" << '\n'
4985 << "/x {lineto stroke} bind def" << '\n'
4986 << "/b {0 0 0 setrgbcolor} def" << '\n'
4987 << "/r {1 0 0 setrgbcolor} def" << '\n';
4988
4989 // calculate colors for level
4990 // coloring; level 0 is black,
4991 // other levels are blue
4992 // ... red
4993 if (eps_flags_base.color_lines_level)
4994 out << "/l { neg " << (max_level) << " add "
4995 << (0.66666 / std::max(1U, (max_level - 1)))
4996 << " mul 1 0.8 sethsbcolor} def" << '\n';
4997
4998 // in 2d, we can also plot cell
4999 // and vertex numbers, but this
5000 // requires a somewhat more
5001 // lengthy preamble. please
5002 // don't ask me what most of
5003 // this means, it is reverse
5004 // engineered from what GNUPLOT
5005 // uses in its output
5006 if ((dim == 2) && (eps_flags_2.write_cell_numbers ||
5007 eps_flags_2.write_vertex_numbers))
5008 {
5009 out
5010 << ("/R {rmoveto} bind def\n"
5011 "/Symbol-Oblique /Symbol findfont [1 0 .167 1 0 0] makefont\n"
5012 "dup length dict begin {1 index /FID eq {pop pop} {def} ifelse} forall\n"
5013 "currentdict end definefont\n"
5014 "/MFshow {{dup dup 0 get findfont exch 1 get scalefont setfont\n"
5015 "[ currentpoint ] exch dup 2 get 0 exch rmoveto dup dup 5 get exch 4 get\n"
5016 "{show} {stringwidth pop 0 rmoveto}ifelse dup 3 get\n"
5017 "{2 get neg 0 exch rmoveto pop} {pop aload pop moveto}ifelse} forall} bind def\n"
5018 "/MFwidth {0 exch {dup 3 get{dup dup 0 get findfont exch 1 get scalefont setfont\n"
5019 "5 get stringwidth pop add}\n"
5020 "{pop} ifelse} forall} bind def\n"
5021 "/MCshow { currentpoint stroke m\n"
5022 "exch dup MFwidth -2 div 3 -1 roll R MFshow } def\n")
5023 << '\n';
5024 }
5025
5026 out << "%%EndProlog" << '\n' << '\n';
5027
5028 // set fine lines
5029 out << eps_flags_base.line_width << " setlinewidth" << '\n';
5030 }
5031
5032 // now write the lines
5033 const Point<2> offset(x_min, y_min);
5034
5035 for (LineList::const_iterator line = line_list.begin();
5036 line != line_list.end();
5037 ++line)
5038 if (eps_flags_base.color_lines_level && (line->level > 0))
5039 out << line->level << " l " << (line->first - offset) * scale << " m "
5040 << (line->second - offset) * scale << " x" << '\n';
5041 else
5042 out << ((line->colorize && eps_flags_base.color_lines_on_user_flag) ?
5043 "r " :
5044 "b ")
5045 << (line->first - offset) * scale << " m "
5046 << (line->second - offset) * scale << " x" << '\n';
5047
5048 // finally write the cell numbers
5049 // in 2d, if that is desired
5050 if ((dim == 2) && (eps_flags_2.write_cell_numbers == true))
5051 {
5052 out << "(Helvetica) findfont 140 scalefont setfont" << '\n';
5053
5054 for (const auto &cell : tria.active_cell_iterators())
5055 {
5056 out << (cell->center()[0] - offset[0]) * scale << ' '
5057 << (cell->center()[1] - offset[1]) * scale << " m" << '\n'
5058 << "[ [(Helvetica) 12.0 0.0 true true (";
5059 if (eps_flags_2.write_cell_number_level)
5060 out << cell;
5061 else
5062 out << cell->index();
5063
5064 out << ")] "
5065 << "] -6 MCshow" << '\n';
5066 }
5067 }
5068
5069 // and the vertex numbers
5070 if ((dim == 2) && (eps_flags_2.write_vertex_numbers == true))
5071 {
5072 out << "(Helvetica) findfont 140 scalefont setfont" << '\n';
5073
5074 // have a list of those
5075 // vertices which we have
5076 // already tracked, to avoid
5077 // doing this multiply
5078 std::set<unsigned int> treated_vertices;
5079 for (const auto &cell : tria.active_cell_iterators())
5080 for (const unsigned int vertex_no : cell->vertex_indices())
5081 if (treated_vertices.find(cell->vertex_index(vertex_no)) ==
5082 treated_vertices.end())
5083 {
5084 treated_vertices.insert(cell->vertex_index(vertex_no));
5085
5086 out << (cell->vertex(vertex_no)[0] - offset[0]) * scale << ' '
5087 << (cell->vertex(vertex_no)[1] - offset[1]) * scale
5088 << " m" << '\n'
5089 << "[ [(Helvetica) 10.0 0.0 true true ("
5090 << cell->vertex_index(vertex_no) << ")] "
5091 << "] -6 MCshow" << '\n';
5092 }
5093 }
5094
5095 out << "showpage" << '\n';
5096
5097 // make sure everything now gets to
5098 // disk
5099 out.flush();
5100
5101 AssertThrow(out.fail() == false, ExcIO());
5102 }
5103 } // namespace
5104} // namespace internal
5105
5106
5107template <int dim, int spacedim>
5108void
5110 std::ostream &out,
5111 const Mapping<dim, spacedim> *mapping) const
5112{
5113 internal::write_eps(tria, out, mapping, eps_flags_2, eps_flags_3);
5114}
5115
5116
5117template <int dim, int spacedim>
5118void
5120 std::ostream &out,
5121 const OutputFormat output_format,
5122 const Mapping<dim, spacedim> *mapping) const
5123{
5124 switch (output_format)
5125 {
5126 case none:
5127 return;
5128
5129 case dx:
5130 write_dx(tria, out);
5131 return;
5132
5133 case ucd:
5134 write_ucd(tria, out);
5135 return;
5136
5137 case gnuplot:
5138 write_gnuplot(tria, out, mapping);
5139 return;
5140
5141 case eps:
5142 write_eps(tria, out, mapping);
5143 return;
5144
5145 case xfig:
5146 write_xfig(tria, out, mapping);
5147 return;
5148
5149 case msh:
5150 write_msh(tria, out);
5151 return;
5152
5153 case svg:
5154 write_svg(tria, out);
5155 return;
5156
5157 case mathgl:
5158 write_mathgl(tria, out);
5159 return;
5160
5161 case vtk:
5162 write_vtk(tria, out);
5163 return;
5164
5165 case vtu:
5166 write_vtu(tria, out);
5167 return;
5168 }
5169
5171}
5172
5173
5174template <int dim, int spacedim>
5175void
5177 std::ostream &out,
5178 const Mapping<dim, spacedim> *mapping) const
5179{
5180 write(tria, out, default_format, mapping);
5181}
5182
5183
5184// explicit instantiations
5185#include "grid_out.inst"
5186
5187
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:3779
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:4100
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:4058
void write_mathgl(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition grid_out.cc:2982
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:3305
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:3938
void write_eps(const Triangulation< dim, spacedim > &tria, std::ostream &out, const Mapping< dim, spacedim > *mapping=nullptr) const
Definition grid_out.cc:5109
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:5119
void write_vtu(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition grid_out.cc:3544
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:3795
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:3899
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:4598
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:3591
@ 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:316
void enter_subsection(const std::string &subsection, const bool create_path_if_needed=true)
long int get_integer(const std::string &entry_string) const
bool get_bool(const std::string &entry_name) const
void declare_entry(const std::string &entry, const std::string &default_value, const Patterns::PatternBase &pattern=Patterns::Anything(), const std::string &documentation="", const bool has_to_be_set=false)
std::string get(const std::string &entry_string) const
double get_double(const std::string &entry_name) const
Definition point.h:111
static 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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Point< 3 > vertices[4]
Point< 2 > second
Definition grid_out.cc:4614
bool colorize
Definition grid_out.cc:4615
Point< 2 > first
Definition grid_out.cc:4613
unsigned int level
Definition grid_out.cc:4616
unsigned int vertex_indices[2]
const unsigned int v0
const unsigned int v1
IteratorRange< active_face_iterator > active_face_iterators() const
IteratorRange< active_cell_iterator > active_cell_iterators() const
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
IteratorRange< cell_iterator > cell_iterators() const
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcInvalidState()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
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)
spacedim & mapping
const Triangulation< dim, spacedim > & tria
if(marked_vertices.size() !=0) for(auto it
for(unsigned int j=best_vertex+1;j< vertices.size();++j) if(vertices_to_use[j])
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr const ReferenceCell Tetrahedron
constexpr const ReferenceCell Quadrilateral
constexpr const ReferenceCell Wedge
constexpr const ReferenceCell Pyramid
constexpr const ReferenceCell Triangle
constexpr const ReferenceCell Hexahedron
constexpr const ReferenceCell Vertex
constexpr const ReferenceCell Line
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:128
std::string encode_base64(const std::vector< unsigned char > &binary_input)
Definition utilities.cc:433
std::string compress(const std::string &input)
Definition utilities.cc:389
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:470
const types::boundary_id invalid_boundary_id
Definition types.h:292
static constexpr double PI
Definition numbers.h:259
const types::boundary_id internal_face_boundary_id
Definition types.h:312
const types::subdomain_id artificial_subdomain_id
Definition types.h:362
static const unsigned int invalid_unsigned_int
Definition types.h:220
const types::manifold_id flat_manifold_id
Definition types.h:325
const InputIterator OutputIterator out
Definition parallel.h:167
const Iterator const std_cxx20::type_identity_t< Iterator > & end
Definition parallel.h:610
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
unsigned int material_id
Definition types.h:167
static void declare_parameters(ParameterHandler &prm)
void parse_parameters(const ParameterHandler &prm)
ReferenceCell reference_cell
Table< 2, float > data
unsigned int n_subdivisions
std::array< Point< spacedim >, GeometryInfo< dim >::vertices_per_cell > vertices
std::array< unsigned int, GeometryInfo< dim >::faces_per_cell > neighbors
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices()
static void declare_parameters(ParameterHandler &param)
Definition grid_out.cc:66
bool write_all_faces
Definition grid_out.h:82
void parse_parameters(ParameterHandler &param)
Definition grid_out.cc:94
bool write_diameter
Definition grid_out.h:71
DX(const bool write_cells=true, const bool write_faces=false, const bool write_diameter=false, const bool write_measure=false, const bool write_all_faces=true)
Definition grid_out.cc:53
static void declare_parameters(ParameterHandler &param)
Definition grid_out.cc:198
unsigned int n_boundary_face_points
Definition grid_out.h:362
void parse_parameters(ParameterHandler &param)
Definition grid_out.cc:231
EpsFlagsBase(const SizeType size_type=width, const unsigned int size=300, const double line_width=0.5, const bool color_lines_on_user_flag=false, const unsigned int n_boundary_face_points=2, const bool color_lines_level=false)
Definition grid_out.cc:182
void parse_parameters(ParameterHandler &param)
Definition grid_out.cc:265
void parse_parameters(ParameterHandler &param)
Definition grid_out.cc:314
void parse_parameters(ParameterHandler &param)
Definition grid_out.cc:357
unsigned int n_extra_curved_line_points
Definition grid_out.h:253
void parse_parameters(ParameterHandler &param)
Definition grid_out.cc:175
Gnuplot(const bool write_cell_number=false, const unsigned int n_extra_curved_line_points=2, const bool curved_inner_cells=false, const bool write_additional_boundary_lines=true)
Definition grid_out.cc:154
static void declare_parameters(ParameterHandler &param)
Definition grid_out.cc:167
bool write_additional_boundary_lines
Definition grid_out.h:272
void parse_parameters(ParameterHandler &param)
Definition grid_out.cc:457
static void declare_parameters(ParameterHandler &param)
Definition grid_out.cc:451
Msh(const bool write_faces=false, const bool write_lines=false)
Definition grid_out.cc:104
void parse_parameters(ParameterHandler &param)
Definition grid_out.cc:118
static void declare_parameters(ParameterHandler &param)
Definition grid_out.cc:110
bool label_level_number
Definition grid_out.h:759
unsigned int height
Definition grid_out.h:661
bool label_level_subdomain_id
Definition grid_out.h:779
bool label_subdomain_id
Definition grid_out.h:774
Background background
Definition grid_out.h:708
unsigned int line_thickness
Definition grid_out.h:672
bool convert_level_number_to_height
Definition grid_out.h:744
float cell_font_scaling
Definition grid_out.h:755
Coloring coloring
Definition grid_out.h:740
Svg(const unsigned int line_thickness=2, const unsigned int boundary_line_thickness=4, const bool margin=true, const Background background=white, const int azimuth_angle=0, const int polar_angle=0, const Coloring coloring=level_number, const bool convert_level_number_to_height=false, const bool label_level_number=false, const bool label_cell_index=false, const bool label_material_id=false, const bool label_subdomain_id=false, const bool draw_colorbar=false, const bool draw_legend=false, const bool label_boundary_id=false)
Definition grid_out.cc:409
@ level_subdomain_id
Convert the level subdomain id into the cell color.
Definition grid_out.h:737
@ subdomain_id
Convert the subdomain id into the cell color.
Definition grid_out.h:735
@ material_id
Convert the material id into the cell color (default)
Definition grid_out.h:731
@ level_number
Convert the level number into the cell color.
Definition grid_out.h:733
unsigned int width
Definition grid_out.h:667
unsigned int boundary_line_thickness
Definition grid_out.h:676
float level_height_factor
Definition grid_out.h:750
Ucd(const bool write_preamble=false, const bool write_faces=false, const bool write_lines=false)
Definition grid_out.cc:125
static void declare_parameters(ParameterHandler &param)
Definition grid_out.cc:136
void parse_parameters(ParameterHandler &param)
Definition grid_out.cc:145
bool output_only_relevant
Definition grid_out.h:893
bool serialize_triangulation
Definition grid_out.h:914
unsigned int n_boundary_face_points
Definition grid_out.h:582
Point< 2 > scaling
Definition grid_out.h:587
@ level_number
Convert the level into the cell color.
Definition grid_out.h:564
@ material_id
Convert the material id into the cell color.
Definition grid_out.h:562
@ level_subdomain_id
Convert the level subdomain id into the cell color.
Definition grid_out.h:568
@ subdomain_id
Convert the global subdomain id into the cell color.
Definition grid_out.h:566
void parse_parameters(ParameterHandler &param)
Definition grid_out.cc:397
static void declare_parameters(ParameterHandler &param)
Definition grid_out.cc:381
Point< 2 > offset
Definition grid_out.h:593
enum GridOutFlags::XFig::Coloring color_by