Reference documentation for deal.II version 9.5.1
\(\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"))
235 size_type = width;
236 else if (param.get("Size by") == std::string("height"))
237 size_type = 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
263
264
265 void
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
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
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
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
Definition tria.cc:12441
void save(Archive &ar, const unsigned int version) const
virtual types::subdomain_id locally_owned_subdomain() const
Definition tria.cc:14520
unsigned int n_active_cells() const
Definition tria.cc:13883
const std::vector< Point< spacedim > > & get_vertices() const
unsigned int n_used_vertices() const
Definition tria.cc:14443
const std::vector< ReferenceCell > & get_reference_cells() const
Definition tria.cc:14821
active_cell_iterator begin_active(const unsigned int level=0) const
Definition tria.cc:13160
const std::vector< bool > & get_used_vertices() const
Definition tria.cc:14452
cell_iterator begin(const unsigned int level=0) const
Definition tria.cc:13139
cell_iterator end() const
Definition tria.cc:13256
#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
Definition tria.cc:13475
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
Definition tria.cc:13382
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition tria.cc:13369
IteratorRange< cell_iterator > cell_iterators() const
Definition tria.cc:13357
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 & get_hypercube()
constexpr const ReferenceCell Vertex
constexpr const ReferenceCell Line
constexpr const ReferenceCell & get_simplex()
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