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