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