deal.II version GIT relicensing-2206-gaa53ff9447 2024-12-02 09:10:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
grid_out.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 1999 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
18#include <deal.II/base/point.h>
21
23
24#include <deal.II/fe/mapping.h>
25
27#include <deal.II/grid/tria.h>
30
32
33#include <boost/archive/binary_oarchive.hpp>
34
35#ifdef DEAL_II_GMSH_WITH_API
36# include <gmsh.h>
37#endif
38
39#include <algorithm>
40#include <cmath>
41#include <cstring>
42#include <ctime>
43#include <fstream>
44#include <iomanip>
45#include <list>
46#include <set>
47
49
50
51namespace GridOutFlags
52{
53 DX::DX(const bool write_cells,
54 const bool write_faces,
55 const bool write_diameter,
56 const bool write_measure,
57 const bool write_all_faces)
58 : write_cells(write_cells)
59 , write_faces(write_faces)
60 , write_diameter(write_diameter)
61 , write_measure(write_measure)
62 , write_all_faces(write_all_faces)
63 {}
64
65 void
67 {
68 param.declare_entry("Write cells",
69 "true",
71 "Write the mesh connectivity as DX grid cells");
72 param.declare_entry("Write faces",
73 "false",
75 "Write faces of cells. These may be boundary faces "
76 "or all faces between mesh cells, according to "
77 "\"Write all faces\"");
78 param.declare_entry("Write diameter",
79 "false",
81 "If cells are written, additionally write their"
82 " diameter as data for visualization");
83 param.declare_entry("Write measure",
84 "false",
86 "Write the volume of each cell as data");
87 param.declare_entry("Write all faces",
88 "true",
90 "Write all faces, not only boundary");
91 }
92
93 void
95 {
96 write_cells = param.get_bool("Write cells");
97 write_faces = param.get_bool("Write faces");
98 write_diameter = param.get_bool("Write diameter");
99 write_measure = param.get_bool("Write measure");
100 write_all_faces = param.get_bool("Write all faces");
101 }
102
103
104 Msh::Msh(const bool write_faces, const bool write_lines)
105 : write_faces(write_faces)
106 , write_lines(write_lines)
107 {}
108
109 void
111 {
112 param.declare_entry("Write faces", "false", Patterns::Bool());
113 param.declare_entry("Write lines", "false", Patterns::Bool());
114 }
115
116
117 void
119 {
120 write_faces = param.get_bool("Write faces");
121 write_lines = param.get_bool("Write lines");
122 }
123
124
125 Ucd::Ucd(const bool write_preamble,
126 const bool write_faces,
127 const bool write_lines)
128 : write_preamble(write_preamble)
129 , write_faces(write_faces)
130 , write_lines(write_lines)
131 {}
132
133
134
135 void
137 {
138 param.declare_entry("Write preamble", "true", Patterns::Bool());
139 param.declare_entry("Write faces", "false", Patterns::Bool());
140 param.declare_entry("Write lines", "false", Patterns::Bool());
141 }
142
143
144 void
146 {
147 write_preamble = param.get_bool("Write preamble");
148 write_faces = param.get_bool("Write faces");
149 write_lines = param.get_bool("Write lines");
150 }
151
152
153
154 Gnuplot::Gnuplot(const bool write_cell_numbers,
155 const unsigned int n_extra_curved_line_points,
156 const bool curved_inner_cells,
157 const bool write_additional_boundary_lines)
158 : write_cell_numbers(write_cell_numbers)
159 , n_extra_curved_line_points(n_extra_curved_line_points)
160 , curved_inner_cells(curved_inner_cells)
161 , write_additional_boundary_lines(write_additional_boundary_lines)
162 {}
163
164
165
166 void
168 {
169 param.declare_entry("Cell number", "false", Patterns::Bool());
170 param.declare_entry("Boundary points", "2", Patterns::Integer());
171 }
172
173
174 void
176 {
177 write_cell_numbers = param.get_bool("Cell number");
178 n_extra_curved_line_points = param.get_integer("Boundary points");
179 }
180
181
183 const unsigned int size,
184 const double line_width,
185 const bool color_lines_on_user_flag,
186 const unsigned int n_boundary_face_points,
187 const bool color_lines_level)
188 : size_type(size_type)
189 , size(size)
190 , line_width(line_width)
191 , color_lines_on_user_flag(color_lines_on_user_flag)
192 , n_boundary_face_points(n_boundary_face_points)
193 , color_lines_level(color_lines_level)
194 {}
195
196
197 void
199 {
200 param.declare_entry("Size by",
201 "width",
202 Patterns::Selection("width|height"),
203 "Depending on this parameter, either the "
204 "width or height "
205 "of the eps is scaled to \"Size\"");
206 param.declare_entry("Size",
207 "300",
209 "Size of the output in points");
210 param.declare_entry("Line width",
211 "0.5",
213 "Width of the lines drawn in points");
214 param.declare_entry("Color by flag",
215 "false",
217 "Draw lines with user flag set in different color");
218 param.declare_entry("Boundary points",
219 "2",
221 "Number of points on boundary edges. "
222 "Increase this beyond 2 to see curved boundaries.");
223 param.declare_entry("Color by level",
224 "false",
226 "Draw different colors according to grid level.");
227 }
228
229
230 void
232 {
233 if (param.get("Size by") == "width")
235 else if (param.get("Size by") == "height")
237 size = param.get_integer("Size");
238 line_width = param.get_double("Line width");
239 color_lines_on_user_flag = param.get_bool("Color by flag");
240 n_boundary_face_points = param.get_integer("Boundary points");
241 color_lines_level = param.get_bool("Color by level");
242 }
243
244
245
246 Eps<1>::Eps(const SizeType size_type,
247 const unsigned int size,
248 const double line_width,
249 const bool color_lines_on_user_flag,
250 const unsigned int n_boundary_face_points)
251 : EpsFlagsBase(size_type,
252 size,
253 line_width,
254 color_lines_on_user_flag,
255 n_boundary_face_points)
256 {}
257
258
259 void
262
263
264 void
269
270
271
272 Eps<2>::Eps(const SizeType size_type,
273 const unsigned int size,
274 const double line_width,
275 const bool color_lines_on_user_flag,
276 const unsigned int n_boundary_face_points,
277 const bool write_cell_numbers,
278 const bool write_cell_number_level,
279 const bool write_vertex_numbers,
280 const bool color_lines_level)
281 : EpsFlagsBase(size_type,
282 size,
283 line_width,
284 color_lines_on_user_flag,
285 n_boundary_face_points,
286 color_lines_level)
287 , write_cell_numbers(write_cell_numbers)
288 , write_cell_number_level(write_cell_number_level)
289 , write_vertex_numbers(write_vertex_numbers)
290 {}
291
292
293 void
295 {
296 param.declare_entry("Cell number",
297 "false",
299 "(2d only) Write cell numbers"
300 " into the centers of cells");
301 param.declare_entry("Level number",
302 "false",
304 "(2d only) if \"Cell number\" is true, write "
305 "numbers in the form level.number");
306 param.declare_entry("Vertex number",
307 "false",
309 "Write numbers for each vertex");
310 }
311
312
313 void
315 {
317 write_cell_numbers = param.get_bool("Cell number");
318 write_cell_number_level = param.get_bool("Level number");
319 write_vertex_numbers = param.get_bool("Vertex number");
320 }
321
322
323
324 Eps<3>::Eps(const SizeType size_type,
325 const unsigned int size,
326 const double line_width,
327 const bool color_lines_on_user_flag,
328 const unsigned int n_boundary_face_points,
329 const double azimut_angle,
330 const double turn_angle)
331 : EpsFlagsBase(size_type,
332 size,
333 line_width,
334 color_lines_on_user_flag,
335 n_boundary_face_points)
336 , azimut_angle(azimut_angle)
337 , turn_angle(turn_angle)
338 {}
339
340
341 void
343 {
344 param.declare_entry("Azimuth",
345 "30",
347 "Azimuth of the view point, that is, the angle "
348 "in the plane from the x-axis.");
349 param.declare_entry("Elevation",
350 "30",
352 "Elevation of the view point above the xy-plane.");
353 }
354
355
356 void
358 {
360 azimut_angle = 90 - param.get_double("Elevation");
361 turn_angle = param.get_double("Azimuth");
362 }
363
364
365
367 : draw_boundary(true)
368 , color_by(material_id)
369 , level_depth(true)
370 , n_boundary_face_points(0)
371 , scaling(1., 1.)
372 , fill_style(20)
373 , line_style(0)
374 , line_thickness(1)
375 , boundary_style(0)
376 , boundary_thickness(3)
377 {}
378
379
380 void
382 {
383 param.declare_entry("Boundary", "true", Patterns::Bool());
384 param.declare_entry("Level color", "false", Patterns::Bool());
385 param.declare_entry("Level depth", "true", Patterns::Bool());
386 // TODO: Unify this number with other output formats
387 param.declare_entry("Boundary points", "0", Patterns::Integer());
388 param.declare_entry("Fill style", "20", Patterns::Integer());
389 param.declare_entry("Line style", "0", Patterns::Integer());
390 param.declare_entry("Line width", "1", Patterns::Integer());
391 param.declare_entry("Boundary style", "0", Patterns::Integer());
392 param.declare_entry("Boundary width", "3", Patterns::Integer());
393 }
394
395
396 void
398 {
399 draw_boundary = param.get_bool("Boundary");
400 level_depth = param.get_bool("Level depth");
401 n_boundary_face_points = param.get_integer("Boundary points");
402 fill_style = param.get_integer("Fill style");
403 line_style = param.get_integer("Line style");
404 line_thickness = param.get_integer("Line width");
405 boundary_style = param.get_integer("Boundary style");
406 boundary_thickness = param.get_integer("Boundary width");
407 }
408
409 Svg::Svg(const unsigned int line_thickness,
410 const unsigned int boundary_line_thickness,
411 bool margin,
412 const Background background,
413 const int azimuth_angle,
414 const int polar_angle,
415 const Coloring coloring,
416 const bool convert_level_number_to_height,
417 const bool label_level_number,
418 const bool label_cell_index,
419 const bool label_material_id,
420 const bool label_subdomain_id,
421 const bool draw_colorbar,
422 const bool draw_legend,
423 const bool label_boundary_id)
424 : height(1000)
425 , width(0)
426 , line_thickness(line_thickness)
427 , boundary_line_thickness(boundary_line_thickness)
428 , margin(margin)
429 , background(background)
430 , azimuth_angle(azimuth_angle)
431 , polar_angle(polar_angle)
432 , coloring(coloring)
433 , convert_level_number_to_height(convert_level_number_to_height)
434 , level_height_factor(0.3f)
435 , cell_font_scaling(1.f)
436 , label_level_number(label_level_number)
437 , label_cell_index(label_cell_index)
438 , label_material_id(label_material_id)
439 , label_subdomain_id(label_subdomain_id)
440 , label_level_subdomain_id(false)
441 , label_boundary_id(label_boundary_id)
442 , draw_colorbar(draw_colorbar)
443 , draw_legend(draw_legend)
444 {}
445
447 : draw_bounding_box(false) // box
448 {}
449
450 void
452 {
453 param.declare_entry("Draw bounding box", "false", Patterns::Bool());
454 }
455
456 void
458 {
459 draw_bounding_box = param.get_bool("Draw bounding box");
460 }
461} // end namespace GridOutFlags
462
463
464
466 : default_format(none)
467{}
468
469
470void
472{
473 dx_flags = flags;
474}
475
476
477
478void
480{
481 msh_flags = flags;
482}
483
484
485void
487{
488 ucd_flags = flags;
489}
490
491
492
493void
495{
496 gnuplot_flags = flags;
497}
498
499
500
501void
503{
504 eps_flags_1 = flags;
505}
506
507
508
509void
511{
512 eps_flags_2 = flags;
513}
514
515
516
517void
519{
520 eps_flags_3 = flags;
521}
522
523
524
525void
527{
528 xfig_flags = flags;
529}
530
531
532void
534{
535 svg_flags = flags;
536}
537
538
539void
541{
542 mathgl_flags = flags;
543}
544
545void
547{
548 vtk_flags = flags;
549}
550
551void
553{
554 vtu_flags = flags;
555}
556
557std::string
559{
560 switch (output_format)
561 {
562 case none:
563 return "";
564 case dx:
565 return ".dx";
566 case gnuplot:
567 return ".gnuplot";
568 case ucd:
569 return ".inp";
570 case eps:
571 return ".eps";
572 case xfig:
573 return ".fig";
574 case msh:
575 return ".msh";
576 case svg:
577 return ".svg";
578 case mathgl:
579 return ".mathgl";
580 case vtk:
581 return ".vtk";
582 case vtu:
583 return ".vtu";
584 default:
586 return "";
587 }
588}
589
590
591
592std::string
597
598
599
601GridOut::parse_output_format(const std::string &format_name)
602{
603 if (format_name == "none" || format_name == "false")
604 return none;
605
606 if (format_name == "dx")
607 return dx;
608
609 if (format_name == "ucd")
610 return ucd;
611
612 if (format_name == "gnuplot")
613 return gnuplot;
614
615 if (format_name == "eps")
616 return eps;
617
618 if (format_name == "xfig")
619 return xfig;
620
621 if (format_name == "msh")
622 return msh;
623
624 if (format_name == "svg")
625 return svg;
626
627 if (format_name == "mathgl")
628 return mathgl;
629
630 if (format_name == "vtk")
631 return vtk;
632
633 if (format_name == "vtu")
634 return vtu;
635
637 // return something weird
638 return OutputFormat(-1);
639}
640
641
642
643std::string
645{
646 return "none|dx|gnuplot|eps|ucd|xfig|msh|svg|mathgl|vtk|vtu";
647}
648
649
650void
696
697
698
699void
701{
702 default_format = parse_output_format(param.get("Format"));
703
704 param.enter_subsection("DX");
706 param.leave_subsection();
707
708 param.enter_subsection("Msh");
710 param.leave_subsection();
711
712 param.enter_subsection("Ucd");
714 param.leave_subsection();
715
716 param.enter_subsection("Gnuplot");
718 param.leave_subsection();
719
720 param.enter_subsection("Eps");
724 param.leave_subsection();
725
726 param.enter_subsection("XFig");
728 param.leave_subsection();
729
730 param.enter_subsection("MathGL");
732 param.leave_subsection();
733
734 param.enter_subsection("Vtk");
736 param.leave_subsection();
737
738 param.enter_subsection("Vtu");
740 param.leave_subsection();
741}
742
743
744
745std::size_t
747{
748 return (sizeof(dx_flags) + sizeof(msh_flags) + sizeof(ucd_flags) +
749 sizeof(gnuplot_flags) + sizeof(eps_flags_1) + sizeof(eps_flags_2) +
750 sizeof(eps_flags_3) + sizeof(xfig_flags) + sizeof(svg_flags) +
751 sizeof(mathgl_flags) + sizeof(vtk_flags) + sizeof(vtu_flags));
752}
753
754
755
756template <>
757void
758GridOut::write_dx(const Triangulation<1> &, std::ostream &) const
759{
761}
762
763template <>
764void
765GridOut::write_dx(const Triangulation<1, 2> &, std::ostream &) const
766{
768}
769
770template <>
771void
772GridOut::write_dx(const Triangulation<1, 3> &, std::ostream &) const
773{
775}
776
777
778
779template <int dim, int spacedim>
780void
782 std::ostream &out) const
783{
784 // TODO:[GK] allow for boundary faces only
786 AssertThrow(out.fail() == false, ExcIO());
787 // Copied and adapted from write_ucd
788 const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
789 const std::vector<bool> &vertex_used = tria.get_used_vertices();
790
791 const unsigned int n_vertices = tria.n_used_vertices();
792
793 // vertices are implicitly numbered from 0 to
794 // n_vertices-1. we have to renumber the
795 // vertices, because otherwise we would end
796 // up with wrong results, if there are unused
797 // vertices
798 std::vector<unsigned int> renumber(vertices.size());
799 // fill this vector with new vertex numbers
800 // ranging from 0 to n_vertices-1
801 unsigned int new_number = 0;
802 for (unsigned int i = 0; i < vertices.size(); ++i)
803 if (vertex_used[i])
804 renumber[i] = new_number++;
805 Assert(new_number == n_vertices, ExcInternalError());
806
807 // write the vertices
808 out << "object \"vertices\" class array type float rank 1 shape " << dim
809 << " items " << n_vertices << " data follows" << '\n';
810
811 for (unsigned int i = 0; i < vertices.size(); ++i)
812 if (vertex_used[i])
813 out << '\t' << vertices[i] << '\n';
814
815 // write cells or faces
816 const bool write_cells = dx_flags.write_cells;
817 const bool write_faces = (dim > 1) ? dx_flags.write_faces : false;
818
819 const unsigned int n_cells = tria.n_active_cells();
820 const unsigned int n_faces =
822
823 const unsigned int n_vertices_per_cell = GeometryInfo<dim>::vertices_per_cell;
824 const unsigned int n_vertices_per_face = GeometryInfo<dim>::vertices_per_face;
825
826 if (write_cells)
827 {
828 out << "object \"cells\" class array type int rank 1 shape "
829 << n_vertices_per_cell << " items " << n_cells << " data follows"
830 << '\n';
831
832 for (const auto &cell : tria.active_cell_iterators())
833 {
834 for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
835 out
836 << '\t'
837 << renumber[cell->vertex_index(GeometryInfo<dim>::dx_to_deal[v])];
838 out << '\n';
839 }
840 out << "attribute \"element type\" string \"";
841 if (dim == 1)
842 out << "lines";
843 if (dim == 2)
844 out << "quads";
845 if (dim == 3)
846 out << "cubes";
847 out << "\"" << '\n'
848 << "attribute \"ref\" string \"positions\"" << '\n'
849 << '\n';
850
851 // Additional cell information
852
853 out << "object \"material\" class array type int rank 0 items " << n_cells
854 << " data follows" << '\n';
855 for (const auto &cell : tria.active_cell_iterators())
856 out << ' ' << cell->material_id();
857 out << '\n' << "attribute \"dep\" string \"connections\"" << '\n' << '\n';
858
859 out << "object \"level\" class array type int rank 0 items " << n_cells
860 << " data follows" << '\n';
861 for (const auto &cell : tria.active_cell_iterators())
862 out << ' ' << cell->level();
863 out << '\n' << "attribute \"dep\" string \"connections\"" << '\n' << '\n';
864
866 {
867 out << "object \"measure\" class array type float rank 0 items "
868 << n_cells << " data follows" << '\n';
869 for (const auto &cell : tria.active_cell_iterators())
870 out << '\t' << cell->measure();
871 out << '\n'
872 << "attribute \"dep\" string \"connections\"" << '\n'
873 << '\n';
874 }
875
877 {
878 out << "object \"diameter\" class array type float rank 0 items "
879 << n_cells << " data follows" << '\n';
880 for (const auto &cell : tria.active_cell_iterators())
881 out << '\t' << cell->diameter();
882 out << '\n'
883 << "attribute \"dep\" string \"connections\"" << '\n'
884 << '\n';
885 }
886 }
887
888 if (write_faces)
889 {
890 out << "object \"faces\" class array type int rank 1 shape "
891 << n_vertices_per_face << " items " << n_faces << " data follows"
892 << '\n';
893
894 for (const auto &cell : tria.active_cell_iterators())
895 {
896 for (const unsigned int f : cell->face_indices())
897 {
899 cell->face(f);
900
901 for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face;
902 ++v)
903 out << '\t'
904 << renumber[face->vertex_index(
906 out << '\n';
907 }
908 }
909 out << "attribute \"element type\" string \"";
910 if (dim == 2)
911 out << "lines";
912 if (dim == 3)
913 out << "quads";
914 out << "\"" << '\n'
915 << "attribute \"ref\" string \"positions\"" << '\n'
916 << '\n';
917
918
919 // Additional face information
920
921 out << "object \"boundary\" class array type int rank 0 items " << n_faces
922 << " data follows" << '\n';
923 for (const auto &cell : tria.active_cell_iterators())
924 {
925 // Little trick to get -1 for the interior
926 for (const unsigned int f : GeometryInfo<dim>::face_indices())
927 {
928 out << ' '
929 << static_cast<std::make_signed_t<types::boundary_id>>(
930 cell->face(f)->boundary_id());
931 }
932 out << '\n';
933 }
934 out << "attribute \"dep\" string \"connections\"" << '\n' << '\n';
935
937 {
938 out << "object \"face measure\" class array type float rank 0 items "
939 << n_faces << " data follows" << '\n';
940 for (const auto &cell : tria.active_cell_iterators())
941 {
942 for (const unsigned int f : GeometryInfo<dim>::face_indices())
943 out << ' ' << cell->face(f)->measure();
944 out << '\n';
945 }
946 out << "attribute \"dep\" string \"connections\"" << '\n' << '\n';
947 }
948
950 {
951 out << "object \"face diameter\" class array type float rank 0 items "
952 << n_faces << " data follows" << '\n';
953 for (const auto &cell : tria.active_cell_iterators())
954 {
955 for (const unsigned int f : GeometryInfo<dim>::face_indices())
956 out << ' ' << cell->face(f)->diameter();
957 out << '\n';
958 }
959 out << "attribute \"dep\" string \"connections\"" << '\n' << '\n';
960 }
961 }
962
963
964 // Write additional face information
965
966 if (write_faces)
967 {
968 }
969 else
970 {
971 }
972
973 // The wrapper
974 out << "object \"deal data\" class field" << '\n'
975 << "component \"positions\" value \"vertices\"" << '\n'
976 << "component \"connections\" value \"cells\"" << '\n';
977
978 if (write_cells)
979 {
980 out << "object \"cell data\" class field" << '\n'
981 << "component \"positions\" value \"vertices\"" << '\n'
982 << "component \"connections\" value \"cells\"" << '\n';
983 out << "component \"material\" value \"material\"" << '\n';
984 out << "component \"level\" value \"level\"" << '\n';
986 out << "component \"measure\" value \"measure\"" << '\n';
988 out << "component \"diameter\" value \"diameter\"" << '\n';
989 }
990
991 if (write_faces)
992 {
993 out << "object \"face data\" class field" << '\n'
994 << "component \"positions\" value \"vertices\"" << '\n'
995 << "component \"connections\" value \"faces\"" << '\n';
996 out << "component \"boundary\" value \"boundary\"" << '\n';
998 out << "component \"measure\" value \"face measure\"" << '\n';
1000 out << "component \"diameter\" value \"face diameter\"" << '\n';
1001 }
1002
1003 out << '\n' << "object \"grid data\" class group" << '\n';
1004 if (write_cells)
1005 out << "member \"cells\" value \"cell data\"" << '\n';
1006 if (write_faces)
1007 out << "member \"faces\" value \"face data\"" << '\n';
1008 out << "end" << '\n';
1009
1010 // make sure everything now gets to
1011 // disk
1012 out.flush();
1013
1014 AssertThrow(out.fail() == false, ExcIO());
1015}
1016
1017
1018
1019template <int dim, int spacedim>
1020void
1022 std::ostream &out) const
1023{
1024 AssertThrow(out.fail() == false, ExcIO());
1025
1026 // get the positions of the
1027 // vertices and whether they are
1028 // used.
1029 const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
1030 const std::vector<bool> &vertex_used = tria.get_used_vertices();
1031
1032 const unsigned int n_vertices = tria.n_used_vertices();
1033
1034 // Write Header
1035 // The file format is:
1036 /*
1037
1038
1039 @f$NOD
1040 number-of-nodes
1041 node-number x-coord y-coord z-coord
1042 ...
1043 @f$ENDNOD
1044 @f$ELM
1045 number-of-elements
1046 elm-number elm-type reg-phys reg-elem number-of-nodes node-number-list
1047 ...
1048 @f$ENDELM
1049 */
1050 out << "@f$NOD" << '\n' << n_vertices << '\n';
1051
1052 // actually write the vertices.
1053 // note that we shall number them
1054 // with first index 1 instead of 0
1055 for (unsigned int i = 0; i < vertices.size(); ++i)
1056 if (vertex_used[i])
1057 {
1058 out << i + 1 // vertex index
1059 << " " << vertices[i];
1060 for (unsigned int d = spacedim + 1; d <= 3; ++d)
1061 out << " 0"; // fill with zeroes
1062 out << '\n';
1063 }
1064
1065 // Write cells preamble
1066 out << "@f$ENDNOD" << '\n'
1067 << "@f$ELM" << '\n'
1068 << tria.n_active_cells() +
1069 ((msh_flags.write_faces ? n_boundary_faces(tria) : 0) +
1071 << '\n';
1072
1073 static constexpr std::array<unsigned int, 8> local_vertex_numbering = {
1074 {0, 1, 5, 4, 2, 3, 7, 6}};
1075
1076 // write cells. Enumerate cells
1077 // consecutively, starting with 1
1078 for (const auto &cell : tria.active_cell_iterators())
1079 {
1080 out << cell->active_cell_index() + 1 << ' '
1081 << cell->reference_cell().gmsh_element_type() << ' '
1082 << cell->material_id() << ' ' << cell->subdomain_id() << ' '
1083 << cell->n_vertices() << ' ';
1084
1085 // Vertex numbering follows UCD conventions.
1086
1087 for (const unsigned int vertex : cell->vertex_indices())
1088 {
1089 if (cell->reference_cell() == ReferenceCells::get_hypercube<dim>())
1090 out << cell->vertex_index(
1091 dim == 3 ? local_vertex_numbering[vertex] :
1093 1
1094 << ' ';
1095 else if (cell->reference_cell() == ReferenceCells::get_simplex<dim>())
1096 out << cell->vertex_index(vertex) + 1 << ' ';
1097 else
1099 }
1100 out << '\n';
1101 }
1102
1103 // write faces and lines with non-zero boundary indicator
1104 unsigned int next_element_index = tria.n_active_cells() + 1;
1106 {
1107 next_element_index = write_msh_faces(tria, next_element_index, out);
1108 }
1110 {
1111 next_element_index = write_msh_lines(tria, next_element_index, out);
1112 }
1113
1114 out << "@f$ENDELM\n";
1115
1116 // make sure everything now gets to
1117 // disk
1118 out.flush();
1119
1120 AssertThrow(out.fail() == false, ExcIO());
1121}
1122
1123
1124template <int dim, int spacedim>
1125void
1127 std::ostream &out) const
1128{
1129 AssertThrow(out.fail() == false, ExcIO());
1130
1131 // get the positions of the
1132 // vertices and whether they are
1133 // used.
1134 const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
1135 const std::vector<bool> &vertex_used = tria.get_used_vertices();
1136
1137 const unsigned int n_vertices = tria.n_used_vertices();
1138
1139 // write preamble
1141 {
1142 // block this to have local
1143 // variables destroyed after
1144 // use
1145 std::time_t time1 = std::time(nullptr);
1146 std::tm *time = std::localtime(&time1);
1147 out
1148 << "# This file was generated by the deal.II library." << '\n'
1149 << "# Date = " << time->tm_year + 1900 << "/" << time->tm_mon + 1
1150 << "/" << time->tm_mday << '\n'
1151 << "# Time = " << time->tm_hour << ":" << std::setw(2) << time->tm_min
1152 << ":" << std::setw(2) << time->tm_sec << '\n'
1153 << "#" << '\n'
1154 << "# For a description of the UCD format see the AVS Developer's guide."
1155 << '\n'
1156 << "#" << '\n';
1157 }
1158
1159 // start with ucd data
1160 out << n_vertices << ' '
1161 << tria.n_active_cells() +
1162 ((ucd_flags.write_faces ? n_boundary_faces(tria) : 0) +
1164 << " 0 0 0" // no data
1165 << '\n';
1166
1167 // actually write the vertices.
1168 // note that we shall number them
1169 // with first index 1 instead of 0
1170 for (unsigned int i = 0; i < vertices.size(); ++i)
1171 if (vertex_used[i])
1172 {
1173 out << i + 1 // vertex index
1174 << " " << vertices[i];
1175 for (unsigned int d = spacedim + 1; d <= 3; ++d)
1176 out << " 0"; // fill with zeroes
1177 out << '\n';
1178 }
1179
1180 // write cells. Enumerate cells
1181 // consecutively, starting with 1
1182 for (const auto &cell : tria.active_cell_iterators())
1183 {
1184 out << cell->active_cell_index() + 1 << ' ' << cell->material_id() << ' ';
1185 switch (dim)
1186 {
1187 case 1:
1188 out << "line ";
1189 break;
1190 case 2:
1191 out << "quad ";
1192 break;
1193 case 3:
1194 out << "hex ";
1195 break;
1196 default:
1198 }
1199
1200 // it follows a list of the
1201 // vertices of each cell. in 1d
1202 // this is simply a list of the
1203 // two vertices, in 2d its counter
1204 // clockwise, as usual in this
1205 // library. in 3d, the same applies
1206 // (special thanks to AVS for
1207 // numbering their vertices in a
1208 // way compatible to deal.II!)
1209 //
1210 // technical reference:
1211 // AVS Developer's Guide, Release 4,
1212 // May, 1992, p. E6
1213 //
1214 // note: vertex numbers are 1-base
1215 for (const unsigned int vertex : GeometryInfo<dim>::vertex_indices())
1216 out << cell->vertex_index(GeometryInfo<dim>::ucd_to_deal[vertex]) + 1
1217 << ' ';
1218 out << '\n';
1219 }
1220
1221 // write faces and lines with non-zero boundary indicator
1222 unsigned int next_element_index = tria.n_active_cells() + 1;
1224 {
1225 next_element_index = write_ucd_faces(tria, next_element_index, out);
1226 }
1228 {
1229 next_element_index = write_ucd_lines(tria, next_element_index, out);
1230 }
1231
1232 // make sure everything now gets to
1233 // disk
1234 out.flush();
1235
1236 AssertThrow(out.fail() == false, ExcIO());
1237}
1238
1239
1240
1241template <int dim, int spacedim>
1242void
1244 std::ostream &,
1245 const Mapping<dim, spacedim> *) const
1246{
1248}
1249
1250
1251// TODO:[GK] Obey parameters
1252template <>
1253void
1255 std::ostream &out,
1256 const Mapping<2> * /*mapping*/) const
1257{
1258 const int dim = 2;
1259 const int spacedim = 2;
1260
1261 const unsigned int nv = GeometryInfo<dim>::vertices_per_cell;
1262
1263 // The following text was copied
1264 // from an existing XFig file.
1265 out << "#FIG 3.2\nLandscape\nCenter\nInches" << std::endl
1266 << "A4\n100.00\nSingle"
1267 << std::endl
1268 // Background is transparent
1269 << "-3" << std::endl
1270 << "# generated by deal.II GridOut class" << std::endl
1271 << "# reduce first number to scale up image" << std::endl
1272 << "1200 2" << std::endl;
1273 // Write custom palette
1274 // grey
1275 unsigned int colno = 32;
1276 out << "0 " << colno++ << " #ff0000" << std::endl;
1277 out << "0 " << colno++ << " #ff8000" << std::endl;
1278 out << "0 " << colno++ << " #ffd000" << std::endl;
1279 out << "0 " << colno++ << " #ffff00" << std::endl;
1280 out << "0 " << colno++ << " #c0ff00" << std::endl;
1281 out << "0 " << colno++ << " #80ff00" << std::endl;
1282 out << "0 " << colno++ << " #00f000" << std::endl;
1283 out << "0 " << colno++ << " #00f0c0" << std::endl;
1284 out << "0 " << colno++ << " #00f0ff" << std::endl;
1285 out << "0 " << colno++ << " #00c0ff" << std::endl;
1286 out << "0 " << colno++ << " #0080ff" << std::endl;
1287 out << "0 " << colno++ << " #0040ff" << std::endl;
1288 out << "0 " << colno++ << " #0000c0" << std::endl;
1289 out << "0 " << colno++ << " #5000ff" << std::endl;
1290 out << "0 " << colno++ << " #8000ff" << std::endl;
1291 out << "0 " << colno++ << " #b000ff" << std::endl;
1292 out << "0 " << colno++ << " #ff00ff" << std::endl;
1293 out << "0 " << colno++ << " #ff80ff" << std::endl;
1294 // grey
1295 for (unsigned int i = 0; i < 8; ++i)
1296 out << "0 " << colno++ << " #" << std::hex << 32 * i + 31 << 32 * i + 31
1297 << 32 * i + 31 << std::dec << std::endl;
1298 // green
1299 for (unsigned int i = 1; i < 16; ++i)
1300 out << "0 " << colno++ << " #00" << std::hex << 16 * i + 15 << std::dec
1301 << "00" << std::endl;
1302 // yellow
1303 for (unsigned int i = 1; i < 16; ++i)
1304 out << "0 " << colno++ << " #" << std::hex << 16 * i + 15 << 16 * i + 15
1305 << std::dec << "00" << std::endl;
1306 // red
1307 for (unsigned int i = 1; i < 16; ++i)
1308 out << "0 " << colno++ << " #" << std::hex << 16 * i + 15 << std::dec
1309 << "0000" << std::endl;
1310 // purple
1311 for (unsigned int i = 1; i < 16; ++i)
1312 out << "0 " << colno++ << " #" << std::hex << 16 * i + 15 << "00"
1313 << 16 * i + 15 << std::dec << std::endl;
1314 // blue
1315 for (unsigned int i = 1; i < 16; ++i)
1316 out << "0 " << colno++ << " #0000" << std::hex << 16 * i + 15 << std::dec
1317 << std::endl;
1318 // cyan
1319 for (unsigned int i = 1; i < 16; ++i)
1320 out << "0 " << colno++ << " #00" << std::hex << 16 * i + 15 << 16 * i + 15
1321 << std::dec << std::endl;
1322
1323 // We write all cells and cells on
1324 // coarser levels are behind cells
1325 // on finer levels. Level 0
1326 // corresponds to a depth of 900,
1327 // each level subtracting 1
1328 for (const auto &cell : tria.cell_iterators())
1329 {
1330 // If depth is not encoded, write finest level only
1331 if (!xfig_flags.level_depth && !cell->is_active())
1332 continue;
1333 // Code for polygon
1334 out << "2 3 " << xfig_flags.line_style << ' '
1336 // with black line
1337 << " 0 ";
1338 // Fill color
1339 switch (xfig_flags.color_by)
1340 {
1341 // TODO[GK]: Simplify after deprecation period is over
1343 out << cell->material_id() + 32;
1344 break;
1346 out << cell->level() + 8;
1347 break;
1349 out << cell->subdomain_id() + 32;
1350 break;
1352 out << cell->level_subdomain_id() + 32;
1353 break;
1354 default:
1356 }
1357
1358 // Depth, unused, fill
1359 out << ' '
1360 << (xfig_flags.level_depth ? (900 - cell->level()) :
1361 (900 + cell->material_id()))
1362 << " 0 " << xfig_flags.fill_style
1363 << " 0.0 "
1364 // some style parameters
1365 << " 0 0 -1 0 0 "
1366 // number of points
1367 << nv + 1 << std::endl;
1368
1369 // For each point, write scaled
1370 // and shifted coordinates
1371 // multiplied by 1200
1372 // (dots/inch)
1373 for (unsigned int k = 0; k <= nv; ++k)
1374 {
1375 const Point<dim> &p =
1376 cell->vertex(GeometryInfo<dim>::ucd_to_deal[k % nv]);
1377 for (unsigned int d = 0; d < static_cast<unsigned int>(dim); ++d)
1378 {
1379 int val = static_cast<int>(1200 * xfig_flags.scaling[d] *
1380 (p[d] - xfig_flags.offset[d]));
1381 out << '\t' << ((d == 0) ? val : -val);
1382 }
1383 out << std::endl;
1384 }
1385 // Now write boundary edges
1386 static const unsigned int face_reorder[4] = {2, 1, 3, 0};
1388 for (const unsigned int f : face_reorder)
1389 {
1391 const types::boundary_id bi = face->boundary_id();
1393 {
1394 // Code for polyline
1395 out << "2 1 "
1396 // with line style and thickness
1397 << xfig_flags.boundary_style << ' '
1398 << xfig_flags.boundary_thickness << ' ' << 1 + bi;
1399 // Fill color
1400 out << " -1 ";
1401 // Depth 100 less than cells
1402 out << (xfig_flags.level_depth ? (800 - cell->level()) :
1403 800 + bi)
1404 // unused, no fill
1405 << " 0 -1 0.0 "
1406 // some style parameters
1407 << " 0 0 -1 0 0 "
1408 // number of points
1410
1411 // For each point, write scaled
1412 // and shifted coordinates
1413 // multiplied by 1200
1414 // (dots/inch)
1415
1416 for (unsigned int k = 0;
1417 k < GeometryInfo<dim>::vertices_per_face;
1418 ++k)
1419 {
1420 const Point<dim> &p = face->vertex(k % nv);
1421 for (unsigned int d = 0; d < static_cast<unsigned int>(dim);
1422 ++d)
1423 {
1424 int val =
1425 static_cast<int>(1200 * xfig_flags.scaling[d] *
1426 (p[d] - xfig_flags.offset[d]));
1427 out << '\t' << ((d == 0) ? val : -val);
1428 }
1429 out << std::endl;
1430 }
1431 }
1432 }
1433 }
1434
1435 // make sure everything now gets to
1436 // disk
1437 out.flush();
1438
1439 AssertThrow(out.fail() == false, ExcIO());
1440}
1441
1442
1443
1444#ifdef DEAL_II_GMSH_WITH_API
1445template <int dim, int spacedim>
1446void
1448 const std::string &filename) const
1449{
1450 // mesh Type renumbering
1451 const std::array<int, 8> dealii_to_gmsh_type = {{15, 1, 2, 3, 4, 7, 6, 5}};
1452
1453 // Vertex renumbering, by dealii type
1454 const std::array<std::vector<unsigned int>, 8> dealii_to_gmsh = {
1455 {{0},
1456 {{0, 1}},
1457 {{0, 1, 2}},
1458 {{0, 1, 3, 2}},
1459 {{0, 1, 2, 3}},
1460 {{0, 1, 3, 2, 4}},
1461 {{0, 1, 2, 3, 4, 5}},
1462 {{0, 1, 3, 2, 4, 5, 7, 6}}}};
1463
1464 // Extract all vertices (nodes in gmsh terminology), and store their three
1465 // dimensional coordinates (regardless of dim).
1466 const auto &vertices = tria.get_vertices();
1467 std::vector<double> coords(3 * vertices.size());
1468 std::vector<std::size_t> nodes(vertices.size());
1469
1470 // Each node has a strictly positive tag. We assign simply its index+1.
1471 std::size_t i = 0;
1472 for (const auto &p : vertices)
1473 {
1474 for (unsigned int d = 0; d < spacedim; ++d)
1475 coords[i * 3 + d] = p[d];
1476 nodes[i] = i + 1;
1477 ++i;
1478 }
1479
1480 // Construct one entity tag per boundary and manifold id pair.
1481 // We need to be smart here, in order to save some disk space. All cells need
1482 // to be written, but only faces and lines that have non default boundary ids
1483 // and/or manifold ids. We collect them into pairs, and for each unique pair,
1484 // we create a gmsh entity where we store the elements. Pre-count all the
1485 // entities, and make sure we know which pair refers to what entity and
1486 // vice-versa.
1487 using IdPair = std::pair<types::material_id, types::manifold_id>;
1488 std::map<IdPair, int> id_pair_to_entity_tag;
1489 std::vector<IdPair> all_pairs;
1490 {
1491 std::set<IdPair> set_of_pairs;
1492 for (const auto &cell : tria.active_cell_iterators())
1493 {
1494 set_of_pairs.insert({cell->material_id(), cell->manifold_id()});
1495 for (const auto &f : cell->face_iterators())
1496 if (f->manifold_id() != numbers::flat_manifold_id ||
1497 (f->boundary_id() != 0 &&
1498 f->boundary_id() != numbers::internal_face_boundary_id))
1499 set_of_pairs.insert({f->boundary_id(), f->manifold_id()});
1500 if (dim > 2)
1501 for (const auto l : cell->line_indices())
1502 {
1503 const auto &f = cell->line(l);
1504 if (f->manifold_id() != numbers::flat_manifold_id ||
1505 (f->boundary_id() != 0 &&
1506 f->boundary_id() != numbers::internal_face_boundary_id))
1507 set_of_pairs.insert({f->boundary_id(), f->manifold_id()});
1508 }
1509 }
1510 all_pairs = {set_of_pairs.begin(), set_of_pairs.end()};
1511
1512 int entity = 1;
1513 for (const auto &p : set_of_pairs)
1514 id_pair_to_entity_tag[p] = entity++;
1515 }
1516
1517 const auto n_entity_tags = id_pair_to_entity_tag.size();
1518
1519 // All elements in the mesh, by entity tag, and by dealii type.
1520 std::vector<std::vector<std::vector<std::size_t>>> element_ids(
1521 n_entity_tags, std::vector<std::vector<std::size_t>>(8));
1522 std::vector<std::vector<std::vector<std::size_t>>> element_nodes(
1523 n_entity_tags, std::vector<std::vector<std::size_t>>(8));
1524
1525 // One element id counter for all dimensions.
1526 std::size_t element_id = 1;
1527
1528 const auto add_element = [&](const auto &element, const int &entity_tag) {
1529 const auto type = element->reference_cell();
1530
1531 Assert(entity_tag > 0, ExcInternalError());
1532 // Add all vertex ids. Make sure we renumber to gmsh, and we add 1 to the
1533 // global index.
1534 for (const auto v : element->vertex_indices())
1535 element_nodes[entity_tag - 1][type].emplace_back(
1536 element->vertex_index(dealii_to_gmsh[type][v]) + 1);
1537
1538 // Save the element id.
1539 element_ids[entity_tag - 1][type].emplace_back(element_id);
1540 ++element_id;
1541 };
1542
1543 // Will create a separate gmsh entity, only if it's a cell, or if the
1544 // boundary and/or the manifold ids are not the default ones.
1545 // In the meanwhile, also store each pair of dimension and entity tag that was
1546 // requested.
1547 std::set<std::pair<int, int>> dim_entity_tag;
1548
1549 auto maybe_add_element =
1550 [&](const auto &element,
1551 const types::boundary_id &boundary_or_material_id) {
1552 const auto struct_dim = element->structure_dimension;
1553 const auto manifold_id = element->manifold_id();
1554
1555 // Exclude default boundary/manifold id or invalid/flag
1556 const bool non_default_boundary_or_material_id =
1557 (boundary_or_material_id != 0 &&
1558 boundary_or_material_id != numbers::internal_face_boundary_id);
1559 const bool non_default_manifold =
1560 manifold_id != numbers::flat_manifold_id;
1561 if (struct_dim == dim || non_default_boundary_or_material_id ||
1562 non_default_manifold)
1563 {
1564 const auto entity_tag =
1565 id_pair_to_entity_tag[{boundary_or_material_id, manifold_id}];
1566 add_element(element, entity_tag);
1567 dim_entity_tag.insert({struct_dim, entity_tag});
1568 }
1569 };
1570
1571 // Loop recursively over all cells, faces, and possibly lines.
1572 for (const auto &cell : tria.active_cell_iterators())
1573 {
1574 maybe_add_element(cell, cell->material_id());
1575 for (const auto &face : cell->face_iterators())
1576 maybe_add_element(face, face->boundary_id());
1577 if (dim > 2)
1578 for (const auto l : cell->line_indices())
1579 maybe_add_element(cell->line(l), cell->line(l)->boundary_id());
1580 }
1581
1582 // Now that we collected everything, plug them into gmsh
1583 gmsh::initialize();
1584 gmsh::option::setNumber("General.Verbosity", 0);
1585 gmsh::model::add("Grid generated in deal.II");
1586 for (const auto &p : dim_entity_tag)
1587 {
1588 gmsh::model::addDiscreteEntity(p.first, p.second);
1589 gmsh::model::mesh::addNodes(p.first, p.second, nodes, coords);
1590 }
1591
1592 for (unsigned int entity_tag = 0; entity_tag < n_entity_tags; ++entity_tag)
1593 for (unsigned int t = 1; t < 8; ++t)
1594 {
1595 const auto all_element_ids = element_ids[entity_tag][t];
1596 const auto all_element_nodes = element_nodes[entity_tag][t];
1597 const auto gmsh_t = dealii_to_gmsh_type[t];
1598 if (all_element_ids.size() > 0)
1599 gmsh::model::mesh::addElementsByType(entity_tag + 1,
1600 gmsh_t,
1601 all_element_ids,
1602 all_element_nodes);
1603 }
1604
1605 // Now for each individual pair of dim and entry, add a physical group, if
1606 // necessary
1607 for (const auto &it : dim_entity_tag)
1608 {
1609 const auto &d = it.first;
1610 const auto &entity_tag = it.second;
1611 const auto &boundary_id = all_pairs[entity_tag - 1].first;
1612 const auto &manifold_id = all_pairs[entity_tag - 1].second;
1613
1614 std::string physical_name;
1615 if (d == dim && boundary_id != 0)
1616 physical_name += "MaterialID:" + Utilities::int_to_string(
1617 static_cast<int>(boundary_id));
1618 else if (d < dim && boundary_id != 0)
1619 physical_name +=
1620 "BoundaryID:" +
1621 (boundary_id == numbers::internal_face_boundary_id ?
1622 "-1" :
1623 Utilities::int_to_string(static_cast<int>(boundary_id)));
1624
1625 std::string sep = physical_name != "" ? ", " : "";
1626 if (manifold_id != numbers::flat_manifold_id)
1627 physical_name +=
1628 sep + "ManifoldID:" +
1629 Utilities::int_to_string(static_cast<int>(manifold_id));
1630 const auto physical_tag =
1631 gmsh::model::addPhysicalGroup(d, {entity_tag}, -1);
1632 if (physical_name != "")
1633 gmsh::model::setPhysicalName(d, physical_tag, physical_name);
1634 }
1635
1636
1637 gmsh::write(filename);
1638 gmsh::clear();
1639 gmsh::finalize();
1640}
1641#endif
1642
1643
1644
1645namespace
1646{
1657 Point<2>
1658 svg_project_point(const Point<3> &point,
1659 const Point<3> &camera_position,
1660 const Tensor<1, 3> &camera_direction,
1661 const Tensor<1, 3> &camera_horizontal,
1662 const float camera_focus)
1663 {
1664 const Tensor<1, 3> camera_vertical =
1665 cross_product_3d(camera_horizontal, camera_direction);
1666
1667 const float phi =
1668 camera_focus / ((point - camera_position) * camera_direction);
1669
1670 const Point<3> projection =
1671 camera_position + phi * (point - camera_position);
1672
1673 return {(projection - camera_position - camera_focus * camera_direction) *
1674 camera_horizontal,
1675 (projection - camera_position - camera_focus * camera_direction) *
1676 camera_vertical};
1677 }
1678} // namespace
1679
1680
1681
1682template <int dim, int spacedim>
1683void
1685 std::ostream & /*out*/) const
1686{
1687 Assert(false,
1688 ExcMessage("Mesh output in SVG format is not implemented for anything "
1689 "other than two-dimensional meshes in two-dimensional "
1690 "space. That's because three-dimensional meshes are best "
1691 "viewed in programs that allow changing the viewpoint, "
1692 "but SVG format does not allow this: It is an inherently "
1693 "2d format, and for three-dimensional meshes would "
1694 "require choosing one, fixed viewpoint."
1695 "\n\n"
1696 "You probably want to output your mesh in a format such "
1697 "as VTK, VTU, or gnuplot."));
1698}
1699
1700
1701void
1702GridOut::write_svg(const Triangulation<2, 2> &tria, std::ostream &out) const
1703{
1704 unsigned int n = 0;
1705
1706 unsigned int min_level, max_level;
1707
1708 // Svg files require an underlying drawing grid. The size of this
1709 // grid is provided in the parameters height and width. Each of them
1710 // may be zero, such that it is computed from the other. Obviously,
1711 // both of them zero does not produce reasonable output.
1712 unsigned int height = svg_flags.height;
1713 unsigned int width = svg_flags.width;
1714 Assert(height != 0 || width != 0,
1715 ExcMessage("You have to set at least one of width and height"));
1716
1717 unsigned int margin_in_percent = 0;
1719 margin_in_percent = 8;
1720
1721 // initial font size for cell labels
1722 unsigned int cell_label_font_size;
1723
1724 // get date and time
1725 // time_t time_stamp;
1726 // tm *now;
1727 // time_stamp = time(0);
1728 // now = localtime(&time_stamp);
1729
1730 float camera_focus;
1731
1732 Point<3> point;
1733 Point<2> projection_decomposition;
1734
1735 float x_max_perspective, x_min_perspective;
1736 float y_max_perspective, y_min_perspective;
1737
1738 float x_dimension_perspective, y_dimension_perspective;
1739
1740
1741 // auxiliary variables for the bounding box and the range of cell levels
1742 double x_min = tria.begin()->vertex(0)[0];
1743 double x_max = x_min;
1744 double y_min = tria.begin()->vertex(0)[1];
1745 double y_max = y_min;
1746
1747 double x_dimension, y_dimension;
1748
1749 min_level = max_level = tria.begin()->level();
1750
1751 // auxiliary set for the materials being used
1752 std::set<unsigned int> materials;
1753
1754 // auxiliary set for the levels being used
1755 std::set<unsigned int> levels;
1756
1757 // auxiliary set for the subdomains being used
1758 std::set<unsigned int> subdomains;
1759
1760 // auxiliary set for the level subdomains being used
1761 std::set<int> level_subdomains;
1762
1763 // We use an active cell iterator to determine the
1764 // bounding box of the given triangulation and check
1765 // the cells for material id, level number, subdomain id
1766 // (, and level subdomain id).
1767 for (const auto &cell : tria.cell_iterators())
1768 {
1769 for (unsigned int vertex_index = 0; vertex_index < cell->n_vertices();
1770 ++vertex_index)
1771 {
1772 if (cell->vertex(vertex_index)[0] < x_min)
1773 x_min = cell->vertex(vertex_index)[0];
1774 if (cell->vertex(vertex_index)[0] > x_max)
1775 x_max = cell->vertex(vertex_index)[0];
1776
1777 if (cell->vertex(vertex_index)[1] < y_min)
1778 y_min = cell->vertex(vertex_index)[1];
1779 if (cell->vertex(vertex_index)[1] > y_max)
1780 y_max = cell->vertex(vertex_index)[1];
1781 }
1782
1783 if (static_cast<unsigned int>(cell->level()) < min_level)
1784 min_level = cell->level();
1785 if (static_cast<unsigned int>(cell->level()) > max_level)
1786 max_level = cell->level();
1787
1788 materials.insert(cell->material_id());
1789 levels.insert(cell->level());
1790 if (cell->is_active())
1791 subdomains.insert(cell->subdomain_id() + 2);
1792 level_subdomains.insert(cell->level_subdomain_id() + 2);
1793 }
1794
1795 x_dimension = x_max - x_min;
1796 y_dimension = y_max - y_min;
1797
1798 // count the materials being used
1799 const unsigned int n_materials = materials.size();
1800
1801 // count the levels being used
1802 const unsigned int n_levels = levels.size();
1803
1804 // count the subdomains being used
1805 const unsigned int n_subdomains = subdomains.size();
1806
1807 // count the level subdomains being used
1808 const unsigned int n_level_subdomains = level_subdomains.size();
1809
1810 switch (svg_flags.coloring)
1811 {
1813 n = n_materials;
1814 break;
1816 n = n_levels;
1817 break;
1819 n = n_subdomains;
1820 break;
1822 n = n_level_subdomains;
1823 break;
1824 default:
1825 break;
1826 }
1827
1828 // set the camera position to top view, targeting at the origin
1829 // vectors and variables for the perspective view
1830 Point<3> camera_position;
1831 camera_position[0] = 0;
1832 camera_position[1] = 0;
1833 camera_position[2] = 2. * std::max(x_dimension, y_dimension);
1834
1835 Tensor<1, 3> camera_direction;
1836 camera_direction[0] = 0;
1837 camera_direction[1] = 0;
1838 camera_direction[2] = -1;
1839
1840 Tensor<1, 3> camera_horizontal;
1841 camera_horizontal[0] = 1;
1842 camera_horizontal[1] = 0;
1843 camera_horizontal[2] = 0;
1844
1845 camera_focus = .5 * std::max(x_dimension, y_dimension);
1846
1847 Point<3> camera_position_temp;
1848 Point<3> camera_direction_temp;
1849 Point<3> camera_horizontal_temp;
1850
1851 const double angle_factor = 3.14159265 / 180.;
1852
1853 // (I) rotate the camera to the chosen polar angle
1854 camera_position_temp[1] =
1855 std::cos(angle_factor * svg_flags.polar_angle) * camera_position[1] -
1856 std::sin(angle_factor * svg_flags.polar_angle) * camera_position[2];
1857 camera_position_temp[2] =
1858 std::sin(angle_factor * svg_flags.polar_angle) * camera_position[1] +
1859 std::cos(angle_factor * svg_flags.polar_angle) * camera_position[2];
1860
1861 camera_direction_temp[1] =
1862 std::cos(angle_factor * svg_flags.polar_angle) * camera_direction[1] -
1863 std::sin(angle_factor * svg_flags.polar_angle) * camera_direction[2];
1864 camera_direction_temp[2] =
1865 std::sin(angle_factor * svg_flags.polar_angle) * camera_direction[1] +
1866 std::cos(angle_factor * svg_flags.polar_angle) * camera_direction[2];
1867
1868 camera_horizontal_temp[1] =
1869 std::cos(angle_factor * svg_flags.polar_angle) * camera_horizontal[1] -
1870 std::sin(angle_factor * svg_flags.polar_angle) * camera_horizontal[2];
1871 camera_horizontal_temp[2] =
1872 std::sin(angle_factor * svg_flags.polar_angle) * camera_horizontal[1] +
1873 std::cos(angle_factor * svg_flags.polar_angle) * camera_horizontal[2];
1874
1875 camera_position[1] = camera_position_temp[1];
1876 camera_position[2] = camera_position_temp[2];
1877
1878 camera_direction[1] = camera_direction_temp[1];
1879 camera_direction[2] = camera_direction_temp[2];
1880
1881 camera_horizontal[1] = camera_horizontal_temp[1];
1882 camera_horizontal[2] = camera_horizontal_temp[2];
1883
1884 // (II) rotate the camera to the chosen azimuth angle
1885 camera_position_temp[0] =
1886 std::cos(angle_factor * svg_flags.azimuth_angle) * camera_position[0] -
1887 std::sin(angle_factor * svg_flags.azimuth_angle) * camera_position[1];
1888 camera_position_temp[1] =
1889 std::sin(angle_factor * svg_flags.azimuth_angle) * camera_position[0] +
1890 std::cos(angle_factor * svg_flags.azimuth_angle) * camera_position[1];
1891
1892 camera_direction_temp[0] =
1893 std::cos(angle_factor * svg_flags.azimuth_angle) * camera_direction[0] -
1894 std::sin(angle_factor * svg_flags.azimuth_angle) * camera_direction[1];
1895 camera_direction_temp[1] =
1896 std::sin(angle_factor * svg_flags.azimuth_angle) * camera_direction[0] +
1897 std::cos(angle_factor * svg_flags.azimuth_angle) * camera_direction[1];
1898
1899 camera_horizontal_temp[0] =
1900 std::cos(angle_factor * svg_flags.azimuth_angle) * camera_horizontal[0] -
1901 std::sin(angle_factor * svg_flags.azimuth_angle) * camera_horizontal[1];
1902 camera_horizontal_temp[1] =
1903 std::sin(angle_factor * svg_flags.azimuth_angle) * camera_horizontal[0] +
1904 std::cos(angle_factor * svg_flags.azimuth_angle) * camera_horizontal[1];
1905
1906 camera_position[0] = camera_position_temp[0];
1907 camera_position[1] = camera_position_temp[1];
1908
1909 camera_direction[0] = camera_direction_temp[0];
1910 camera_direction[1] = camera_direction_temp[1];
1911
1912 camera_horizontal[0] = camera_horizontal_temp[0];
1913 camera_horizontal[1] = camera_horizontal_temp[1];
1914
1915 // translate the camera to the given triangulation
1916 camera_position[0] = x_min + .5 * x_dimension;
1917 camera_position[1] = y_min + .5 * y_dimension;
1918
1919 camera_position[0] += 2. * std::max(x_dimension, y_dimension) *
1920 std::sin(angle_factor * svg_flags.polar_angle) *
1921 std::sin(angle_factor * svg_flags.azimuth_angle);
1922 camera_position[1] -= 2. * std::max(x_dimension, y_dimension) *
1923 std::sin(angle_factor * svg_flags.polar_angle) *
1924 std::cos(angle_factor * svg_flags.azimuth_angle);
1925
1926
1927 // determine the bounding box of the given triangulation on the projection
1928 // plane of the camera viewing system
1929 point[0] = tria.begin()->vertex(0)[0];
1930 point[1] = tria.begin()->vertex(0)[1];
1931 point[2] = 0;
1932
1933 float min_level_min_vertex_distance = 0;
1934
1936 {
1937 point[2] = svg_flags.level_height_factor *
1938 (static_cast<float>(tria.begin()->level()) /
1939 static_cast<float>(n_levels)) *
1940 std::max(x_dimension, y_dimension);
1941 }
1942
1943 projection_decomposition = svg_project_point(
1944 point, camera_position, camera_direction, camera_horizontal, camera_focus);
1945
1946 x_max_perspective = projection_decomposition[0];
1947 x_min_perspective = projection_decomposition[0];
1948
1949 y_max_perspective = projection_decomposition[1];
1950 y_min_perspective = projection_decomposition[1];
1951
1952 for (const auto &cell : tria.cell_iterators())
1953 {
1954 point[0] = cell->vertex(0)[0];
1955 point[1] = cell->vertex(0)[1];
1956 point[2] = 0;
1957
1959 {
1960 point[2] =
1962 (static_cast<float>(cell->level()) / static_cast<float>(n_levels)) *
1963 std::max(x_dimension, y_dimension);
1964 }
1965
1966 projection_decomposition = svg_project_point(point,
1967 camera_position,
1968 camera_direction,
1969 camera_horizontal,
1970 camera_focus);
1971
1972 if (x_max_perspective < projection_decomposition[0])
1973 x_max_perspective = projection_decomposition[0];
1974 if (x_min_perspective > projection_decomposition[0])
1975 x_min_perspective = projection_decomposition[0];
1976
1977 if (y_max_perspective < projection_decomposition[1])
1978 y_max_perspective = projection_decomposition[1];
1979 if (y_min_perspective > projection_decomposition[1])
1980 y_min_perspective = projection_decomposition[1];
1981
1982 point[0] = cell->vertex(1)[0];
1983 point[1] = cell->vertex(1)[1];
1984
1985 projection_decomposition = svg_project_point(point,
1986 camera_position,
1987 camera_direction,
1988 camera_horizontal,
1989 camera_focus);
1990
1991 if (x_max_perspective < projection_decomposition[0])
1992 x_max_perspective = projection_decomposition[0];
1993 if (x_min_perspective > projection_decomposition[0])
1994 x_min_perspective = projection_decomposition[0];
1995
1996 if (y_max_perspective < projection_decomposition[1])
1997 y_max_perspective = projection_decomposition[1];
1998 if (y_min_perspective > projection_decomposition[1])
1999 y_min_perspective = projection_decomposition[1];
2000
2001 point[0] = cell->vertex(2)[0];
2002 point[1] = cell->vertex(2)[1];
2003
2004 projection_decomposition = svg_project_point(point,
2005 camera_position,
2006 camera_direction,
2007 camera_horizontal,
2008 camera_focus);
2009
2010 if (x_max_perspective < projection_decomposition[0])
2011 x_max_perspective = projection_decomposition[0];
2012 if (x_min_perspective > projection_decomposition[0])
2013 x_min_perspective = projection_decomposition[0];
2014
2015 if (y_max_perspective < projection_decomposition[1])
2016 y_max_perspective = projection_decomposition[1];
2017 if (y_min_perspective > projection_decomposition[1])
2018 y_min_perspective = projection_decomposition[1];
2019
2020 if (cell->n_vertices() == 4) // in case of quadrilateral
2021 {
2022 point[0] = cell->vertex(3)[0];
2023 point[1] = cell->vertex(3)[1];
2024
2025 projection_decomposition = svg_project_point(point,
2026 camera_position,
2027 camera_direction,
2028 camera_horizontal,
2029 camera_focus);
2030
2031 if (x_max_perspective < projection_decomposition[0])
2032 x_max_perspective = projection_decomposition[0];
2033 if (x_min_perspective > projection_decomposition[0])
2034 x_min_perspective = projection_decomposition[0];
2035
2036 if (y_max_perspective < projection_decomposition[1])
2037 y_max_perspective = projection_decomposition[1];
2038 if (y_min_perspective > projection_decomposition[1])
2039 y_min_perspective = projection_decomposition[1];
2040 }
2041
2042 if (static_cast<unsigned int>(cell->level()) == min_level)
2043 min_level_min_vertex_distance = cell->minimum_vertex_distance();
2044 }
2045
2046 x_dimension_perspective = x_max_perspective - x_min_perspective;
2047 y_dimension_perspective = y_max_perspective - y_min_perspective;
2048
2049 // create the svg file with an internal style sheet
2050 if (width == 0)
2051 width = static_cast<unsigned int>(
2052 .5 + height * (x_dimension_perspective / y_dimension_perspective));
2053 else if (height == 0)
2054 height = static_cast<unsigned int>(
2055 .5 + width * (y_dimension_perspective / x_dimension_perspective));
2056 unsigned int additional_width = 0;
2057 // font size for date, time, legend, and colorbar
2058 unsigned int font_size =
2059 static_cast<unsigned int>(.5 + (height / 100.) * 1.75);
2060 cell_label_font_size = static_cast<unsigned int>(
2061 .5 + (height * .15 * svg_flags.cell_font_scaling *
2062 min_level_min_vertex_distance / std::min(x_dimension, y_dimension)));
2063
2064 if (svg_flags.draw_legend &&
2068 {
2069 additional_width = static_cast<unsigned int>(
2070 .5 + height * .4); // additional width for legend
2071 }
2072 else if (svg_flags.draw_colorbar && (svg_flags.coloring != 0u))
2073 {
2074 additional_width = static_cast<unsigned int>(
2075 .5 + height * .175); // additional width for colorbar
2076 }
2077
2078 // out << "<!-- deal.ii GridOut " << now->tm_mday << '/' << now->tm_mon + 1 <<
2079 // '/' << now->tm_year + 1900
2080 // << ' ' << now->tm_hour << ':';
2081 //
2082 // if (now->tm_min < 10) out << '0';
2083 //
2084 // out << now->tm_min << " -->" << '\n';
2085
2086 // basic svg header
2087 out << "<svg width=\"" << width + additional_width << "\" height=\"" << height
2088 << "\" xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\">" << '\n'
2089 << '\n';
2090
2091
2093 {
2094 out
2095 << " <linearGradient id=\"background_gradient\" gradientUnits=\"userSpaceOnUse\" x1=\"0\" y1=\"0\" x2=\"0\" y2=\""
2096 << height << "\">" << '\n'
2097 << " <stop offset=\"0\" style=\"stop-color:white\"/>" << '\n'
2098 << " <stop offset=\"1\" style=\"stop-color:lightsteelblue\"/>" << '\n'
2099 << " </linearGradient>" << '\n';
2100 }
2101
2102 out << '\n';
2103
2104 // header for the internal style sheet
2105 out << "<!-- internal style sheet -->" << '\n'
2106 << "<style type=\"text/css\"><![CDATA[" << '\n';
2107
2108 // set the background of the output graphic
2110 out << " rect.background{fill:url(#background_gradient)}" << '\n';
2112 out << " rect.background{fill:white}" << '\n';
2113 else
2114 out << " rect.background{fill:none}" << '\n';
2115
2116 // basic svg graphic element styles
2117 out << " rect{fill:none; stroke:rgb(25,25,25); stroke-width:"
2118 << svg_flags.line_thickness << '}' << '\n'
2119 << " text{font-family:Helvetica; text-anchor:middle; fill:rgb(25,25,25)}"
2120 << '\n'
2121 << " line{stroke:rgb(25,25,25); stroke-width:"
2122 << svg_flags.boundary_line_thickness << '}' << '\n'
2123 << " path{fill:none; stroke:rgb(25,25,25); stroke-width:"
2124 << svg_flags.line_thickness << '}' << '\n'
2125 << " circle{fill:white; stroke:black; stroke-width:2}" << '\n'
2126 << '\n';
2127
2128 // polygon styles with respect to the chosen cell coloring
2129 if (svg_flags.coloring != 0u)
2130 {
2131 unsigned int labeling_index = 0;
2132 auto materials_it = materials.begin();
2133 auto levels_it = levels.begin();
2134 auto subdomains_it = subdomains.begin();
2135 auto level_subdomains_it = level_subdomains.begin();
2136
2137 for (unsigned int index = 0; index < n; ++index)
2138 {
2139 double h;
2140
2141 if (n != 1)
2142 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;
3182 tria.save_user_flags_line(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;
3224 tria.save_user_flags_line(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_mpi_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 const 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 const std::vector<double> dummy_weights(n_points, 1. / n_points);
4237 const 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 auto offset =
4292 cell->reference_cell(),
4293 face_no,
4294 cell->combined_face_orientation(face_no),
4295 n_points);
4296 for (unsigned int i = 0; i < n_points; ++i)
4297 line_points.push_back(
4299 cell, q_projector.point(offset + i)));
4300 internal::remove_colinear_points(line_points);
4301
4302 for (const Point<spacedim> &point : line_points)
4303 out << point << ' ' << cell->level() << ' '
4304 << cell->material_id() << '\n';
4305
4306 out << '\n' << '\n';
4307 }
4308 else
4309 {
4310 // if, however, the face is not at the boundary and we
4311 // don't want to curve anything, then draw it as usual
4312 out << face->vertex(0) << ' ' << cell->level() << ' '
4313 << cell->material_id() << '\n'
4314 << face->vertex(1) << ' ' << cell->level() << ' '
4315 << cell->material_id() << '\n'
4316 << '\n'
4317 << '\n';
4318 }
4319 }
4320 }
4321 }
4322
4323 // make sure everything now gets to disk
4324 out.flush();
4325
4326 AssertThrow(out.fail() == false, ExcIO());
4327 }
4328
4329
4330
4331 template <int spacedim>
4332 void
4333 write_gnuplot(const ::Triangulation<3, spacedim> &tria,
4334 std::ostream &out,
4335 const Mapping<3, spacedim> *mapping,
4336 const GridOutFlags::Gnuplot &gnuplot_flags)
4337 {
4338 AssertThrow(out.fail() == false, ExcIO());
4339
4340 const int dim = 3;
4341
4342 const unsigned int n_additional_points =
4343 gnuplot_flags.n_extra_curved_line_points;
4344 const unsigned int n_points = 2 + n_additional_points;
4345
4346 // If we need to plot curved lines then generate a quadrature formula to
4347 // place points via the mapping
4348 std::unique_ptr<Quadrature<dim>> q_projector;
4349 std::vector<Point<1>> boundary_points;
4350 if (mapping != nullptr)
4351 {
4352 boundary_points.resize(n_points);
4353 boundary_points[0][0] = 0;
4354 boundary_points[n_points - 1][0] = 1;
4355 for (unsigned int i = 1; i < n_points - 1; ++i)
4356 boundary_points[i][0] = 1. * i / (n_points - 1);
4357
4358 const std::vector<double> dummy_weights(n_points, 1. / n_points);
4359 const Quadrature<1> quadrature1d(boundary_points, dummy_weights);
4360
4361 // tensor product of points, only one copy
4362 const QIterated<dim - 1> quadrature(quadrature1d, 1);
4363 q_projector = std::make_unique<Quadrature<dim>>(
4365 ReferenceCells::get_hypercube<dim>(), quadrature));
4366 }
4367
4368 for (const auto &cell : tria.active_cell_iterators())
4369 {
4370 if (gnuplot_flags.write_cell_numbers)
4371 out << "# cell " << cell << '\n';
4372
4373 if (mapping == nullptr || n_points == 2 ||
4374 (!cell->has_boundary_lines() &&
4375 !gnuplot_flags.curved_inner_cells))
4376 {
4377 if (cell->reference_cell() == ReferenceCells::Hexahedron)
4378 {
4379 // front face
4380 out << cell->vertex(0) << ' ' << cell->level() << ' '
4381 << cell->material_id() << '\n'
4382 << cell->vertex(1) << ' ' << cell->level() << ' '
4383 << cell->material_id() << '\n'
4384 << cell->vertex(5) << ' ' << cell->level() << ' '
4385 << cell->material_id() << '\n'
4386 << cell->vertex(4) << ' ' << cell->level() << ' '
4387 << cell->material_id() << '\n'
4388 << cell->vertex(0) << ' ' << cell->level() << ' '
4389 << cell->material_id() << '\n'
4390 << '\n';
4391 // back face
4392 out << cell->vertex(2) << ' ' << cell->level() << ' '
4393 << cell->material_id() << '\n'
4394 << cell->vertex(3) << ' ' << cell->level() << ' '
4395 << cell->material_id() << '\n'
4396 << cell->vertex(7) << ' ' << cell->level() << ' '
4397 << cell->material_id() << '\n'
4398 << cell->vertex(6) << ' ' << cell->level() << ' '
4399 << cell->material_id() << '\n'
4400 << cell->vertex(2) << ' ' << cell->level() << ' '
4401 << cell->material_id() << '\n'
4402 << '\n';
4403
4404 // now for the four connecting lines
4405 out << cell->vertex(0) << ' ' << cell->level() << ' '
4406 << cell->material_id() << '\n'
4407 << cell->vertex(2) << ' ' << cell->level() << ' '
4408 << cell->material_id() << '\n'
4409 << '\n';
4410 out << cell->vertex(1) << ' ' << cell->level() << ' '
4411 << cell->material_id() << '\n'
4412 << cell->vertex(3) << ' ' << cell->level() << ' '
4413 << cell->material_id() << '\n'
4414 << '\n';
4415 out << cell->vertex(5) << ' ' << cell->level() << ' '
4416 << cell->material_id() << '\n'
4417 << cell->vertex(7) << ' ' << cell->level() << ' '
4418 << cell->material_id() << '\n'
4419 << '\n';
4420 out << cell->vertex(4) << ' ' << cell->level() << ' '
4421 << cell->material_id() << '\n'
4422 << cell->vertex(6) << ' ' << cell->level() << ' '
4423 << cell->material_id() << '\n'
4424 << '\n';
4425 }
4426 else if (cell->reference_cell() == ReferenceCells::Tetrahedron)
4427 {
4428 // Draw the tetrahedron as a two collections of lines.
4429 for (const unsigned int v : {0, 1, 2, 0, 3, 2})
4430 {
4431 out << cell->vertex(v) << ' ' << cell->level() << ' '
4432 << cell->material_id() << '\n';
4433 }
4434 out << '\n'; // end of first line
4435
4436 for (const unsigned int v : {3, 1})
4437 {
4438 out << cell->vertex(v) << ' ' << cell->level() << ' '
4439 << cell->material_id() << '\n';
4440 }
4441 out << '\n'; // end of second line
4442 }
4443 else if (cell->reference_cell() == ReferenceCells::Wedge)
4444 {
4445 // Draw the wedge as a collection of three
4446 // lines. The first one wraps around the base,
4447 // goes up to the top, and wraps around that. The
4448 // second and third are just individual lines
4449 // going from base to top.
4450 for (const unsigned int v : {0, 1, 2, 0, 3, 4, 5, 3})
4451 {
4452 out << cell->vertex(v) << ' ' << cell->level() << ' '
4453 << cell->material_id() << '\n';
4454 }
4455 out << '\n'; // end of first line
4456
4457 for (const unsigned int v : {1, 4})
4458 {
4459 out << cell->vertex(v) << ' ' << cell->level() << ' '
4460 << cell->material_id() << '\n';
4461 }
4462 out << '\n'; // end of second line
4463
4464 for (const unsigned int v : {2, 5})
4465 {
4466 out << cell->vertex(v) << ' ' << cell->level() << ' '
4467 << cell->material_id() << '\n';
4468 }
4469 out << '\n'; // end of third line
4470 }
4471 else if (cell->reference_cell() == ReferenceCells::Pyramid)
4472 {
4473 // Draw the pyramid as a collections of two lines.
4474 for (const unsigned int v : {0, 1, 3, 2, 0, 4, 1})
4475 {
4476 out << cell->vertex(v) << ' ' << cell->level() << ' '
4477 << cell->material_id() << '\n';
4478 }
4479 out << '\n'; // end of first line
4480
4481 for (const unsigned int v : {2, 4, 3})
4482 {
4483 out << cell->vertex(v) << ' ' << cell->level() << ' '
4484 << cell->material_id() << '\n';
4485 }
4486 out << '\n'; // end of second line
4487 }
4488 else
4490 }
4491 else // need to handle curved boundaries
4492 {
4493 Assert(cell->reference_cell() == ReferenceCells::Hexahedron,
4495 for (const unsigned int face_no :
4496 GeometryInfo<dim>::face_indices())
4497 {
4498 const typename ::Triangulation<dim,
4499 spacedim>::face_iterator
4500 face = cell->face(face_no);
4501
4502 if (face->at_boundary() &&
4503 gnuplot_flags.write_additional_boundary_lines)
4504 {
4505 const auto offset =
4507 cell->reference_cell(),
4508 face_no,
4509 cell->combined_face_orientation(face_no),
4510 n_points * n_points);
4511 for (unsigned int i = 0; i < n_points - 1; ++i)
4512 for (unsigned int j = 0; j < n_points - 1; ++j)
4513 {
4514 const Point<spacedim> p0 =
4516 cell,
4517 q_projector->point(offset + i * n_points + j));
4518 out << p0 << ' ' << cell->level() << ' '
4519 << cell->material_id() << '\n';
4520 out << (mapping->transform_unit_to_real_cell(
4521 cell,
4522 q_projector->point(
4523 offset + (i + 1) * n_points + j)))
4524 << ' ' << cell->level() << ' '
4525 << cell->material_id() << '\n';
4526 out << (mapping->transform_unit_to_real_cell(
4527 cell,
4528 q_projector->point(
4529 offset + (i + 1) * n_points + j + 1)))
4530 << ' ' << cell->level() << ' '
4531 << cell->material_id() << '\n';
4532 out << (mapping->transform_unit_to_real_cell(
4533 cell,
4534 q_projector->point(offset + i * n_points +
4535 j + 1)))
4536 << ' ' << cell->level() << ' '
4537 << cell->material_id() << '\n';
4538 // and the first point again
4539 out << p0 << ' ' << cell->level() << ' '
4540 << cell->material_id() << '\n';
4541 out << '\n' << '\n';
4542 }
4543 }
4544 else
4545 {
4546 for (unsigned int l = 0;
4547 l < GeometryInfo<dim>::lines_per_face;
4548 ++l)
4549 {
4550 const typename ::Triangulation<dim, spacedim>::
4551 line_iterator line = face->line(l);
4552
4553 const Point<spacedim> &v0 = line->vertex(0),
4554 &v1 = line->vertex(1);
4555 if (line->at_boundary() ||
4556 gnuplot_flags.curved_inner_cells)
4557 {
4558 // Save the points on each face to a vector and
4559 // then try to remove colinear points that won't
4560 // show up in the generated plot.
4561 std::vector<Point<spacedim>> line_points;
4562 // transform_real_to_unit_cell could be replaced
4563 // by using QProjector<dim>::project_to_line
4564 // which is not yet implemented
4565 const Point<spacedim>
4566 u0 = mapping->transform_real_to_unit_cell(cell,
4567 v0),
4568 u1 = mapping->transform_real_to_unit_cell(cell,
4569 v1);
4570 for (unsigned int i = 0; i < n_points; ++i)
4571 line_points.push_back(
4573 cell,
4574 (1 - boundary_points[i][0]) * u0 +
4575 boundary_points[i][0] * u1));
4576 internal::remove_colinear_points(line_points);
4577 for (const Point<spacedim> &point : line_points)
4578 out << point << ' ' << cell->level() << ' '
4579 << static_cast<unsigned int>(
4580 cell->material_id())
4581 << '\n';
4582 }
4583 else
4584 out << v0 << ' ' << cell->level() << ' '
4585 << cell->material_id() << '\n'
4586 << v1 << ' ' << cell->level() << ' '
4587 << cell->material_id() << '\n';
4588
4589 out << '\n' << '\n';
4590 }
4591 }
4592 }
4593 }
4594 }
4595
4596 // make sure everything now gets to disk
4597 out.flush();
4598
4599 AssertThrow(out.fail() == false, ExcIO());
4600 }
4601 } // namespace
4602} // namespace internal
4603
4604
4605
4606template <int dim, int spacedim>
4607void
4609 std::ostream &out,
4610 const Mapping<dim, spacedim> *mapping) const
4611{
4612 internal::write_gnuplot(tria, out, mapping, gnuplot_flags);
4613}
4614
4615
4616
4617namespace internal
4618{
4619 namespace
4620 {
4621 struct LineEntry
4622 {
4626 unsigned int level;
4627 LineEntry(const Point<2> &f,
4628 const Point<2> &s,
4629 const bool c,
4630 const unsigned int l)
4631 : first(f)
4632 , second(s)
4633 , colorize(c)
4634 , level(l)
4635 {}
4636 };
4637
4638
4639 void
4640 write_eps(const ::Triangulation<1> &,
4641 std::ostream &,
4642 const Mapping<1> *,
4643 const GridOutFlags::Eps<2> &,
4644 const GridOutFlags::Eps<3> &)
4645 {
4647 }
4648
4649 void
4650 write_eps(const ::Triangulation<1, 2> &,
4651 std::ostream &,
4652 const Mapping<1, 2> *,
4653 const GridOutFlags::Eps<2> &,
4654 const GridOutFlags::Eps<3> &)
4655 {
4657 }
4658
4659 void
4660 write_eps(const ::Triangulation<1, 3> &,
4661 std::ostream &,
4662 const Mapping<1, 3> *,
4663 const GridOutFlags::Eps<2> &,
4664 const GridOutFlags::Eps<3> &)
4665 {
4667 }
4668
4669 void
4670 write_eps(const ::Triangulation<2, 3> &,
4671 std::ostream &,
4672 const Mapping<2, 3> *,
4673 const GridOutFlags::Eps<2> &,
4674 const GridOutFlags::Eps<3> &)
4675 {
4677 }
4678
4679
4680
4681 template <int dim, int spacedim>
4682 void
4683 write_eps(const ::Triangulation<dim, spacedim> &tria,
4684 std::ostream &out,
4685 const Mapping<dim, spacedim> *mapping,
4686 const GridOutFlags::Eps<2> &eps_flags_2,
4687 const GridOutFlags::Eps<3> &eps_flags_3)
4688 {
4689 using LineList = std::list<LineEntry>;
4690
4691 // We should never get here in 1d since this function is overloaded for
4692 // all dim == 1 cases.
4693 Assert(dim == 2 || dim == 3, ExcInternalError());
4694
4695 // Copy, with an object slice, something containing the flags common to
4696 // all dimensions in order to avoid the recurring distinctions between
4697 // the different eps_flags present.
4698 const GridOutFlags::EpsFlagsBase eps_flags_base =
4699 dim == 2 ?
4700 static_cast<const GridOutFlags::EpsFlagsBase &>(eps_flags_2) :
4701 static_cast<const GridOutFlags::EpsFlagsBase &>(eps_flags_3);
4702
4703 AssertThrow(out.fail() == false, ExcIO());
4704 const unsigned int n_points = eps_flags_base.n_boundary_face_points;
4705
4706 // make up a list of lines by which
4707 // we will construct the triangulation
4708 //
4709 // this part unfortunately is a bit
4710 // dimension dependent, so we have to
4711 // treat every dimension different.
4712 // however, by directly producing
4713 // the lines to be printed, i.e. their
4714 // 2d images, we can later do the
4715 // actual output dimension independent
4716 // again
4717 LineList line_list;
4718
4719 switch (dim)
4720 {
4721 case 1:
4722 {
4724 break;
4725 }
4726
4727 case 2:
4728 {
4729 for (const auto &cell : tria.active_cell_iterators())
4730 for (const unsigned int line_no : cell->line_indices())
4731 {
4732 typename ::Triangulation<dim, spacedim>::line_iterator
4733 line = cell->line(line_no);
4734
4735 // first treat all
4736 // interior lines and
4737 // make up a list of
4738 // them. if curved
4739 // lines shall not be
4740 // supported (i.e. no
4741 // mapping is
4742 // provided), then also
4743 // treat all other
4744 // lines
4745 if (!line->has_children() &&
4746 (mapping == nullptr || !line->at_boundary()))
4747 // one would expect
4748 // make_pair(line->vertex(0),
4749 // line->vertex(1))
4750 // here, but that is
4751 // not dimension
4752 // independent, since
4753 // vertex(i) is
4754 // Point<dim>, but we
4755 // want a Point<2>.
4756 // in fact, whenever
4757 // we're here, the
4758 // vertex is a
4759 // Point<dim>, but
4760 // the compiler does
4761 // not know
4762 // this. hopefully,
4763 // the compiler will
4764 // optimize away this
4765 // little kludge
4766 line_list.emplace_back(
4767 Point<2>(line->vertex(0)[0], line->vertex(0)[1]),
4768 Point<2>(line->vertex(1)[0], line->vertex(1)[1]),
4769 line->user_flag_set(),
4770 cell->level());
4771 }
4772
4773 // next if we are to treat
4774 // curved boundaries
4775 // specially, then add lines
4776 // to the list consisting of
4777 // pieces of the boundary
4778 // lines
4779 if (mapping != nullptr)
4780 {
4781 // to do so, first
4782 // generate a sequence of
4783 // points on a face and
4784 // project them onto the
4785 // faces of a unit cell
4786 std::vector<Point<dim - 1>> boundary_points(n_points);
4787
4788 for (unsigned int i = 0; i < n_points; ++i)
4789 boundary_points[i][0] = 1. * (i + 1) / (n_points + 1);
4790
4791 const Quadrature<dim - 1> quadrature(boundary_points);
4792 const Quadrature<dim> q_projector(
4794 ReferenceCells::get_hypercube<dim>(), quadrature));
4795
4796 // next loop over all
4797 // boundary faces and
4798 // generate the info from
4799 // them
4800 for (const auto &cell : tria.active_cell_iterators())
4801 for (const unsigned int face_no :
4802 GeometryInfo<dim>::face_indices())
4803 {
4804 const typename ::Triangulation<dim, spacedim>::
4805 face_iterator face = cell->face(face_no);
4806
4807 if (face->at_boundary())
4808 {
4809 Point<dim> p0_dim(face->vertex(0));
4810 Point<2> p0(p0_dim[0], p0_dim[1]);
4811
4812 // loop over
4813 // all pieces
4814 // of the line
4815 // and generate
4816 // line-lets
4817 const auto offset =
4819 cell->reference_cell(),
4820 face_no,
4821 cell->combined_face_orientation(face_no),
4822 n_points);
4823 for (unsigned int i = 0; i < n_points; ++i)
4824 {
4825 const Point<dim> p1_dim(
4827 cell, q_projector.point(offset + i)));
4828 const Point<2> p1(p1_dim[0], p1_dim[1]);
4829
4830 line_list.emplace_back(p0,
4831 p1,
4832 face->user_flag_set(),
4833 cell->level());
4834 p0 = p1;
4835 }
4836
4837 // generate last piece
4838 const Point<dim> p1_dim(face->vertex(1));
4839 const Point<2> p1(p1_dim[0], p1_dim[1]);
4840 line_list.emplace_back(p0,
4841 p1,
4842 face->user_flag_set(),
4843 cell->level());
4844 }
4845 }
4846 }
4847
4848 break;
4849 }
4850
4851 case 3:
4852 {
4853 // curved boundary output
4854 // presently not supported
4855 Assert(mapping == nullptr, ExcNotImplemented());
4856
4857 // loop over all lines and compute their
4858 // projection on the plane perpendicular
4859 // to the direction of sight
4860
4861 // direction of view equals the unit
4862 // vector of the position of the
4863 // spectator to the origin.
4864 //
4865 // we chose here the viewpoint as in
4866 // gnuplot as default.
4867 //
4868 // TODO:[WB] Fix a potential problem with viewing angles in 3d Eps
4869 // GridOut
4870 // note: the following might be wrong
4871 // if one of the base vectors below
4872 // is in direction of the viewer, but
4873 // I am too tired at present to fix
4874 // this
4875 const double pi = numbers::PI;
4876 const double z_angle = eps_flags_3.azimut_angle;
4877 const double turn_angle = eps_flags_3.turn_angle;
4878 const Point<dim> view_direction(
4879 -std::sin(z_angle * 2. * pi / 360.) *
4880 std::sin(turn_angle * 2. * pi / 360.),
4881 +std::sin(z_angle * 2. * pi / 360.) *
4882 std::cos(turn_angle * 2. * pi / 360.),
4883 -std::cos(z_angle * 2. * pi / 360.));
4884
4885 // decide about the two unit vectors
4886 // in this plane. we chose the first one
4887 // to be the projection of the z-axis
4888 // to this plane
4889 const Tensor<1, dim> vector1 =
4890 Point<dim>(0, 0, 1) -
4891 ((Point<dim>(0, 0, 1) * view_direction) * view_direction);
4892 const Tensor<1, dim> unit_vector1 = vector1 / vector1.norm();
4893
4894 // now the third vector is fixed. we
4895 // chose the projection of a more or
4896 // less arbitrary vector to the plane
4897 // perpendicular to the first one
4898 const Tensor<1, dim> vector2 =
4899 (Point<dim>(1, 0, 0) -
4900 ((Point<dim>(1, 0, 0) * view_direction) * view_direction) -
4901 ((Point<dim>(1, 0, 0) * unit_vector1) * unit_vector1));
4902 const Tensor<1, dim> unit_vector2 = vector2 / vector2.norm();
4903
4904
4905 for (const auto &cell : tria.active_cell_iterators())
4906 for (const unsigned int line_no : cell->line_indices())
4907 {
4908 typename ::Triangulation<dim, spacedim>::line_iterator
4909 line = cell->line(line_no);
4910 line_list.emplace_back(
4911 Point<2>(line->vertex(0) * unit_vector2,
4912 line->vertex(0) * unit_vector1),
4913 Point<2>(line->vertex(1) * unit_vector2,
4914 line->vertex(1) * unit_vector1),
4915 line->user_flag_set(),
4916 cell->level());
4917 }
4918
4919 break;
4920 }
4921
4922 default:
4924 }
4925
4926
4927
4928 // find out minimum and maximum x and
4929 // y coordinates to compute offsets
4930 // and scaling factors
4931 double x_min = tria.begin_active()->vertex(0)[0];
4932 double x_max = x_min;
4933 double y_min = tria.begin_active()->vertex(0)[1];
4934 double y_max = y_min;
4935 unsigned int max_level = line_list.begin()->level;
4936
4937 for (LineList::const_iterator line = line_list.begin();
4938 line != line_list.end();
4939 ++line)
4940 {
4941 x_min = std::min(x_min, line->first[0]);
4942 x_min = std::min(x_min, line->second[0]);
4943
4944 x_max = std::max(x_max, line->first[0]);
4945 x_max = std::max(x_max, line->second[0]);
4946
4947 y_min = std::min(y_min, line->first[1]);
4948 y_min = std::min(y_min, line->second[1]);
4949
4950 y_max = std::max(y_max, line->first[1]);
4951 y_max = std::max(y_max, line->second[1]);
4952
4953 max_level = std::max(max_level, line->level);
4954 }
4955
4956 // scale in x-direction such that
4957 // in the output 0 <= x <= 300.
4958 // don't scale in y-direction to
4959 // preserve the shape of the
4960 // triangulation
4961 const double scale =
4962 (eps_flags_base.size /
4963 (eps_flags_base.size_type == GridOutFlags::EpsFlagsBase::width ?
4964 x_max - x_min :
4965 y_min - y_max));
4966
4967
4968 // now write preamble
4969 {
4970 // block this to have local
4971 // variables destroyed after
4972 // use
4973 std::time_t time1 = std::time(nullptr);
4974 std::tm *time = std::localtime(&time1);
4975 out << "%!PS-Adobe-2.0 EPSF-1.2" << '\n'
4976 << "%%Title: deal.II Output" << '\n'
4977 << "%%Creator: the deal.II library" << '\n'
4978 << "%%Creation Date: " << time->tm_year + 1900 << "/"
4979 << time->tm_mon + 1 << "/" << time->tm_mday << " - "
4980 << time->tm_hour << ":" << std::setw(2) << time->tm_min << ":"
4981 << std::setw(2) << time->tm_sec << '\n'
4982 << "%%BoundingBox: "
4983 // lower left corner
4984 << "0 0 "
4985 // upper right corner
4986 << static_cast<unsigned int>(
4987 std::floor(((x_max - x_min) * scale) + 1))
4988 << ' '
4989 << static_cast<unsigned int>(
4990 std::floor(((y_max - y_min) * scale) + 1))
4991 << '\n';
4992
4993 // define some abbreviations to keep
4994 // the output small:
4995 // m=move turtle to
4996 // x=execute line stroke
4997 // b=black pen
4998 // r=red pen
4999 out << "/m {moveto} bind def" << '\n'
5000 << "/x {lineto stroke} bind def" << '\n'
5001 << "/b {0 0 0 setrgbcolor} def" << '\n'
5002 << "/r {1 0 0 setrgbcolor} def" << '\n';
5003
5004 // calculate colors for level
5005 // coloring; level 0 is black,
5006 // other levels are blue
5007 // ... red
5008 if (eps_flags_base.color_lines_level)
5009 out << "/l { neg " << (max_level) << " add "
5010 << (0.66666 / std::max(1U, (max_level - 1)))
5011 << " mul 1 0.8 sethsbcolor} def" << '\n';
5012
5013 // in 2d, we can also plot cell
5014 // and vertex numbers, but this
5015 // requires a somewhat more
5016 // lengthy preamble. please
5017 // don't ask me what most of
5018 // this means, it is reverse
5019 // engineered from what GNUPLOT
5020 // uses in its output
5021 if ((dim == 2) && (eps_flags_2.write_cell_numbers ||
5022 eps_flags_2.write_vertex_numbers))
5023 {
5024 out
5025 << ("/R {rmoveto} bind def\n"
5026 "/Symbol-Oblique /Symbol findfont [1 0 .167 1 0 0] makefont\n"
5027 "dup length dict begin {1 index /FID eq {pop pop} {def} ifelse} forall\n"
5028 "currentdict end definefont\n"
5029 "/MFshow {{dup dup 0 get findfont exch 1 get scalefont setfont\n"
5030 "[ currentpoint ] exch dup 2 get 0 exch rmoveto dup dup 5 get exch 4 get\n"
5031 "{show} {stringwidth pop 0 rmoveto}ifelse dup 3 get\n"
5032 "{2 get neg 0 exch rmoveto pop} {pop aload pop moveto}ifelse} forall} bind def\n"
5033 "/MFwidth {0 exch {dup 3 get{dup dup 0 get findfont exch 1 get scalefont setfont\n"
5034 "5 get stringwidth pop add}\n"
5035 "{pop} ifelse} forall} bind def\n"
5036 "/MCshow { currentpoint stroke m\n"
5037 "exch dup MFwidth -2 div 3 -1 roll R MFshow } def\n")
5038 << '\n';
5039 }
5040
5041 out << "%%EndProlog" << '\n' << '\n';
5042
5043 // set fine lines
5044 out << eps_flags_base.line_width << " setlinewidth" << '\n';
5045 }
5046
5047 // now write the lines
5048 const Point<2> offset(x_min, y_min);
5049
5050 for (LineList::const_iterator line = line_list.begin();
5051 line != line_list.end();
5052 ++line)
5053 if (eps_flags_base.color_lines_level && (line->level > 0))
5054 out << line->level << " l " << (line->first - offset) * scale << " m "
5055 << (line->second - offset) * scale << " x" << '\n';
5056 else
5057 out << ((line->colorize && eps_flags_base.color_lines_on_user_flag) ?
5058 "r " :
5059 "b ")
5060 << (line->first - offset) * scale << " m "
5061 << (line->second - offset) * scale << " x" << '\n';
5062
5063 // finally write the cell numbers
5064 // in 2d, if that is desired
5065 if ((dim == 2) && (eps_flags_2.write_cell_numbers == true))
5066 {
5067 out << "(Helvetica) findfont 140 scalefont setfont" << '\n';
5068
5069 for (const auto &cell : tria.active_cell_iterators())
5070 {
5071 out << (cell->center()[0] - offset[0]) * scale << ' '
5072 << (cell->center()[1] - offset[1]) * scale << " m" << '\n'
5073 << "[ [(Helvetica) 12.0 0.0 true true (";
5074 if (eps_flags_2.write_cell_number_level)
5075 out << cell;
5076 else
5077 out << cell->index();
5078
5079 out << ")] "
5080 << "] -6 MCshow" << '\n';
5081 }
5082 }
5083
5084 // and the vertex numbers
5085 if ((dim == 2) && (eps_flags_2.write_vertex_numbers == true))
5086 {
5087 out << "(Helvetica) findfont 140 scalefont setfont" << '\n';
5088
5089 // have a list of those
5090 // vertices which we have
5091 // already tracked, to avoid
5092 // doing this multiply
5093 std::set<unsigned int> treated_vertices;
5094 for (const auto &cell : tria.active_cell_iterators())
5095 for (const unsigned int vertex_no : cell->vertex_indices())
5096 if (treated_vertices.find(cell->vertex_index(vertex_no)) ==
5097 treated_vertices.end())
5098 {
5099 treated_vertices.insert(cell->vertex_index(vertex_no));
5100
5101 out << (cell->vertex(vertex_no)[0] - offset[0]) * scale << ' '
5102 << (cell->vertex(vertex_no)[1] - offset[1]) * scale
5103 << " m" << '\n'
5104 << "[ [(Helvetica) 10.0 0.0 true true ("
5105 << cell->vertex_index(vertex_no) << ")] "
5106 << "] -6 MCshow" << '\n';
5107 }
5108 }
5109
5110 out << "showpage" << '\n';
5111
5112 // make sure everything now gets to
5113 // disk
5114 out.flush();
5115
5116 AssertThrow(out.fail() == false, ExcIO());
5117 }
5118 } // namespace
5119} // namespace internal
5120
5121
5122template <int dim, int spacedim>
5123void
5125 std::ostream &out,
5126 const Mapping<dim, spacedim> *mapping) const
5127{
5128 internal::write_eps(tria, out, mapping, eps_flags_2, eps_flags_3);
5129}
5130
5131
5132template <int dim, int spacedim>
5133void
5135 std::ostream &out,
5136 const OutputFormat output_format,
5137 const Mapping<dim, spacedim> *mapping) const
5138{
5139 switch (output_format)
5140 {
5141 case none:
5142 return;
5143
5144 case dx:
5145 write_dx(tria, out);
5146 return;
5147
5148 case ucd:
5149 write_ucd(tria, out);
5150 return;
5151
5152 case gnuplot:
5153 write_gnuplot(tria, out, mapping);
5154 return;
5155
5156 case eps:
5157 write_eps(tria, out, mapping);
5158 return;
5159
5160 case xfig:
5161 write_xfig(tria, out, mapping);
5162 return;
5163
5164 case msh:
5165 write_msh(tria, out);
5166 return;
5167
5168 case svg:
5169 write_svg(tria, out);
5170 return;
5171
5172 case mathgl:
5173 write_mathgl(tria, out);
5174 return;
5175
5176 case vtk:
5177 write_vtk(tria, out);
5178 return;
5179
5180 case vtu:
5181 write_vtu(tria, out);
5182 return;
5183 }
5184
5186}
5187
5188
5189template <int dim, int spacedim>
5190void
5192 std::ostream &out,
5193 const Mapping<dim, spacedim> *mapping) const
5194{
5195 write(tria, out, default_format, mapping);
5196}
5197
5198
5199// explicit instantiations
5200#include "grid_out.inst"
5201
5202
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:5124
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:5134
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:4608
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:318
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const =0
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const =0
void enter_subsection(const std::string &subsection, const bool create_path_if_needed=true)
long int get_integer(const std::string &entry_string) const
bool get_bool(const std::string &entry_name) const
void declare_entry(const std::string &entry, const std::string &default_value, const Patterns::PatternBase &pattern=Patterns::Anything(), const std::string &documentation="", const bool has_to_be_set=false)
std::string get(const std::string &entry_string) const
double get_double(const std::string &entry_name) const
Definition point.h:111
static DataSetDescriptor face(const ReferenceCell &reference_cell, const unsigned int face_no, const bool face_orientation, const bool face_flip, const bool face_rotation, const unsigned int n_quadrature_points)
static Quadrature< dim > project_to_all_faces(const ReferenceCell &reference_cell, const hp::QCollection< dim - 1 > &quadrature)
const Point< dim > & point(const unsigned int i) const
unsigned int gmsh_element_type() const
numbers::NumberTraits< Number >::real_type norm() const
void save_user_flags_line(std::ostream &out) const
void save(Archive &ar, const unsigned int version) const
cell_iterator begin(const unsigned int level=0) const
virtual types::subdomain_id locally_owned_subdomain() const
unsigned int n_active_cells() const
const std::vector< Point< spacedim > > & get_vertices() const
unsigned int n_used_vertices() const
const std::vector< ReferenceCell > & get_reference_cells() const
cell_iterator end() const
const std::vector< bool > & get_used_vertices() const
active_cell_iterator begin_active(const unsigned int level=0) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
Point< 2 > second
Definition grid_out.cc:4624
bool colorize
Definition grid_out.cc:4625
Point< 2 > first
Definition grid_out.cc:4623
unsigned int level
Definition grid_out.cc:4626
unsigned int vertex_indices[2]
const unsigned int v0
const unsigned int v1
IteratorRange< active_face_iterator > active_face_iterators() const
IteratorRange< active_cell_iterator > active_cell_iterators() const
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
IteratorRange< cell_iterator > cell_iterators() const
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcInvalidState()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
std::size_t size
Definition mpi.cc:734
const unsigned int n_procs
Definition mpi.cc:924
void write_eps(const std::vector< Patch< 2, spacedim > > &patches, const std::vector< std::string > &data_names, const std::vector< std::tuple< unsigned int, unsigned int, std::string, DataComponentInterpretation::DataComponentInterpretation > > &nonscalar_data_ranges, const EpsFlags &flags, std::ostream &out)
void write_vtu_header(std::ostream &out, const VtkFlags &flags)
void write_vtu(const std::vector< Patch< dim, spacedim > > &patches, const std::vector< std::string > &data_names, const std::vector< std::tuple< unsigned int, unsigned int, std::string, DataComponentInterpretation::DataComponentInterpretation > > &nonscalar_data_ranges, const VtkFlags &flags, std::ostream &out)
void write_vtu_main(const std::vector< Patch< dim, spacedim > > &patches, const std::vector< std::string > &data_names, const std::vector< std::tuple< unsigned int, unsigned int, std::string, DataComponentInterpretation::DataComponentInterpretation > > &nonscalar_data_ranges, const VtkFlags &flags, std::ostream &out)
void write_vtu_footer(std::ostream &out)
void write_gnuplot(const std::vector< Patch< dim, spacedim > > &patches, const std::vector< std::string > &data_names, const std::vector< std::tuple< unsigned int, unsigned int, std::string, DataComponentInterpretation::DataComponentInterpretation > > &nonscalar_data_ranges, const GnuplotFlags &flags, std::ostream &out)
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr const ReferenceCell Tetrahedron
constexpr const ReferenceCell Quadrilateral
constexpr const ReferenceCell Wedge
constexpr const ReferenceCell Pyramid
constexpr const ReferenceCell Triangle
constexpr const ReferenceCell Hexahedron
constexpr const ReferenceCell Vertex
constexpr const ReferenceCell Line
VectorType::value_type * end(VectorType &V)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:92
std::string encode_base64(const std::vector< unsigned char > &binary_input)
Definition utilities.cc:433
std::string compress(const std::string &input)
Definition utilities.cc:389
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:470
const types::boundary_id invalid_boundary_id
Definition types.h:292
static constexpr double PI
Definition numbers.h:254
const types::boundary_id internal_face_boundary_id
Definition types.h:312
const types::subdomain_id artificial_subdomain_id
Definition types.h:362
static const unsigned int invalid_unsigned_int
Definition types.h:220
const types::manifold_id flat_manifold_id
Definition types.h:325
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
unsigned int material_id
Definition types.h:167
static void declare_parameters(ParameterHandler &prm)
void parse_parameters(const ParameterHandler &prm)
ReferenceCell reference_cell
Table< 2, float > data
unsigned int n_subdivisions
std::array< Point< spacedim >, GeometryInfo< dim >::vertices_per_cell > vertices
std::array< unsigned int, GeometryInfo< dim >::faces_per_cell > neighbors
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices()
static void declare_parameters(ParameterHandler &param)
Definition grid_out.cc:66
bool write_all_faces
Definition grid_out.h:82
void parse_parameters(ParameterHandler &param)
Definition grid_out.cc:94
bool write_diameter
Definition grid_out.h:71
DX(const bool write_cells=true, const bool write_faces=false, const bool write_diameter=false, const bool write_measure=false, const bool write_all_faces=true)
Definition grid_out.cc:53
static void declare_parameters(ParameterHandler &param)
Definition grid_out.cc:198
unsigned int n_boundary_face_points
Definition grid_out.h:362
void parse_parameters(ParameterHandler &param)
Definition grid_out.cc:231
EpsFlagsBase(const SizeType size_type=width, const unsigned int size=300, const double line_width=0.5, const bool color_lines_on_user_flag=false, const unsigned int n_boundary_face_points=2, const bool color_lines_level=false)
Definition grid_out.cc:182
void parse_parameters(ParameterHandler &param)
Definition grid_out.cc:265
void parse_parameters(ParameterHandler &param)
Definition grid_out.cc:314
void parse_parameters(ParameterHandler &param)
Definition grid_out.cc:357
unsigned int n_extra_curved_line_points
Definition grid_out.h:253
void parse_parameters(ParameterHandler &param)
Definition grid_out.cc:175
Gnuplot(const bool write_cell_number=false, const unsigned int n_extra_curved_line_points=2, const bool curved_inner_cells=false, const bool write_additional_boundary_lines=true)
Definition grid_out.cc:154
static void declare_parameters(ParameterHandler &param)
Definition grid_out.cc:167
bool write_additional_boundary_lines
Definition grid_out.h:272
void parse_parameters(ParameterHandler &param)
Definition grid_out.cc:457
static void declare_parameters(ParameterHandler &param)
Definition grid_out.cc:451
Msh(const bool write_faces=false, const bool write_lines=false)
Definition grid_out.cc:104
void parse_parameters(ParameterHandler &param)
Definition grid_out.cc:118
static void declare_parameters(ParameterHandler &param)
Definition grid_out.cc:110
bool label_level_number
Definition grid_out.h:759
unsigned int height
Definition grid_out.h:661
bool label_level_subdomain_id
Definition grid_out.h:779
bool label_subdomain_id
Definition grid_out.h:774
Background background
Definition grid_out.h:708
unsigned int line_thickness
Definition grid_out.h:672
bool convert_level_number_to_height
Definition grid_out.h:744
float cell_font_scaling
Definition grid_out.h:755
Coloring coloring
Definition grid_out.h:740
Svg(const unsigned int line_thickness=2, const unsigned int boundary_line_thickness=4, const bool margin=true, const Background background=white, const int azimuth_angle=0, const int polar_angle=0, const Coloring coloring=level_number, const bool convert_level_number_to_height=false, const bool label_level_number=false, const bool label_cell_index=false, const bool label_material_id=false, const bool label_subdomain_id=false, const bool draw_colorbar=false, const bool draw_legend=false, const bool label_boundary_id=false)
Definition grid_out.cc:409
@ level_subdomain_id
Convert the level subdomain id into the cell color.
Definition grid_out.h:737
@ subdomain_id
Convert the subdomain id into the cell color.
Definition grid_out.h:735
@ material_id
Convert the material id into the cell color (default)
Definition grid_out.h:731
@ level_number
Convert the level number into the cell color.
Definition grid_out.h:733
unsigned int width
Definition grid_out.h:667
unsigned int boundary_line_thickness
Definition grid_out.h:676
float level_height_factor
Definition grid_out.h:750