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