34 #include <boost/archive/binary_oarchive.hpp>
36 #ifdef DEAL_II_GMSH_WITH_API
55 const bool write_faces,
56 const bool write_diameter,
57 const bool write_measure,
58 const bool write_all_faces)
60 , write_faces(write_faces)
61 , write_diameter(write_diameter)
62 , write_measure(write_measure)
63 , write_all_faces(write_all_faces)
72 "Write the mesh connectivity as DX grid cells");
76 "Write faces of cells. These may be boundary faces "
77 "or all faces between mesh cells, according to "
78 "\"Write all faces\"");
82 "If cells are written, additionally write their"
83 " diameter as data for visualization");
87 "Write the volume of each cell as data");
91 "Write all faces, not only boundary");
105 Msh::Msh(
const bool write_faces,
const bool write_lines)
106 : write_faces(write_faces)
107 , write_lines(write_lines)
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)
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)
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)
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)
204 "Depending on this parameter, either the "
206 "of the eps is scaled to \"Size\"");
210 "Size of the output in points");
214 "Width of the lines drawn in points");
218 "Draw lines with user flag set in different color");
222 "Number of points on boundary edges. "
223 "Increase this beyond 2 to see curved boundaries.");
227 "Draw different colors according to grid level.");
234 if (param.
get(
"Size by") == std::string(
"width"))
236 else if (param.
get(
"Size by") == std::string(
"height"))
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)
255 color_lines_on_user_flag,
256 n_boundary_face_points)
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)
285 color_lines_on_user_flag,
286 n_boundary_face_points,
288 , write_cell_numbers(write_cell_numbers)
289 , write_cell_number_level(write_cell_number_level)
290 , write_vertex_numbers(write_vertex_numbers)
300 "(2D only) Write cell numbers"
301 " into the centers of cells");
305 "(2D only) if \"Cell number\" is true, write "
306 "numbers in the form level.number");
310 "Write numbers for each vertex");
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");
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)
335 color_lines_on_user_flag,
336 n_boundary_face_points)
337 , azimut_angle(azimut_angle)
338 , turn_angle(turn_angle)
348 "Azimuth of the viw point, that is, the angle "
349 "in the plane from the x-axis.");
353 "Elevation of the view point above the xy-plane.");
361 azimut_angle = 90 - param.
get_double(
"Elevation");
368 : draw_boundary(true)
371 , n_boundary_face_points(0)
377 , boundary_thickness(3)
411 const unsigned int boundary_line_thickness,
414 const int azimuth_angle,
415 const int polar_angle,
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)
427 , line_thickness(line_thickness)
428 , boundary_line_thickness(boundary_line_thickness)
430 , background(background)
431 , azimuth_angle(azimuth_angle)
432 , polar_angle(polar_angle)
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)
448 : draw_bounding_box(false)
561 switch (output_format)
604 if (format_name ==
"none" || format_name ==
"false")
607 if (format_name ==
"dx")
610 if (format_name ==
"ucd")
613 if (format_name ==
"gnuplot")
616 if (format_name ==
"eps")
619 if (format_name ==
"xfig")
622 if (format_name ==
"msh")
625 if (format_name ==
"svg")
628 if (format_name ==
"mathgl")
631 if (format_name ==
"vtk")
634 if (format_name ==
"vtu")
647 return "none|dx|gnuplot|eps|ucd|xfig|msh|svg|mathgl|vtk|vtu";
780 template <
int dim,
int spacedim>
783 std::ostream & out)
const
799 std::vector<unsigned int> renumber(
vertices.size());
802 unsigned int new_number = 0;
803 for (
unsigned int i = 0; i <
vertices.size(); ++i)
805 renumber[i] = new_number++;
809 out <<
"object \"vertices\" class array type float rank 1 shape " << dim
810 <<
" items " << n_vertices <<
" data follows" <<
'\n';
812 for (
unsigned int i = 0; i <
vertices.size(); ++i)
821 const unsigned int n_faces =
829 out <<
"object \"cells\" class array type int rank 1 shape "
830 << n_vertices_per_cell <<
" items " <<
n_cells <<
" data follows"
841 out <<
"attribute \"element type\" string \"";
849 <<
"attribute \"ref\" string \"positions\"" <<
'\n'
854 out <<
"object \"material\" class array type int rank 0 items " <<
n_cells
855 <<
" data follows" <<
'\n';
857 out <<
' ' << cell->material_id();
858 out <<
'\n' <<
"attribute \"dep\" string \"connections\"" <<
'\n' <<
'\n';
860 out <<
"object \"level\" class array type int rank 0 items " <<
n_cells
861 <<
" data follows" <<
'\n';
863 out <<
' ' << cell->level();
864 out <<
'\n' <<
"attribute \"dep\" string \"connections\"" <<
'\n' <<
'\n';
868 out <<
"object \"measure\" class array type float rank 0 items "
869 <<
n_cells <<
" data follows" <<
'\n';
871 out <<
'\t' << cell->measure();
873 <<
"attribute \"dep\" string \"connections\"" <<
'\n'
879 out <<
"object \"diameter\" class array type float rank 0 items "
880 <<
n_cells <<
" data follows" <<
'\n';
882 out <<
'\t' << cell->diameter();
884 <<
"attribute \"dep\" string \"connections\"" <<
'\n'
891 out <<
"object \"faces\" class array type int rank 1 shape "
892 << n_vertices_per_face <<
" items " << n_faces <<
" data follows"
897 for (
const unsigned int f : cell->face_indices())
902 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face;
905 << renumber[face->vertex_index(
910 out <<
"attribute \"element type\" string \"";
916 <<
"attribute \"ref\" string \"positions\"" <<
'\n'
922 out <<
"object \"boundary\" class array type int rank 0 items " << n_faces
923 <<
" data follows" <<
'\n';
930 <<
static_cast<std::make_signed<types::boundary_id>::type
>(
931 cell->face(f)->boundary_id());
935 out <<
"attribute \"dep\" string \"connections\"" <<
'\n' <<
'\n';
939 out <<
"object \"face measure\" class array type float rank 0 items "
940 << n_faces <<
" data follows" <<
'\n';
944 out <<
' ' << cell->face(f)->measure();
947 out <<
"attribute \"dep\" string \"connections\"" <<
'\n' <<
'\n';
952 out <<
"object \"face diameter\" class array type float rank 0 items "
953 << n_faces <<
" data follows" <<
'\n';
957 out <<
' ' << cell->face(f)->diameter();
960 out <<
"attribute \"dep\" string \"connections\"" <<
'\n' <<
'\n';
973 out <<
"object \"deal data\" class field" <<
'\n'
974 <<
"component \"positions\" value \"vertices\"" <<
'\n'
975 <<
"component \"connections\" value \"cells\"" <<
'\n';
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';
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';
1002 out <<
'\n' <<
"object \"grid data\" class group" <<
'\n';
1004 out <<
"member \"cells\" value \"cell data\"" <<
'\n';
1006 out <<
"member \"faces\" value \"face data\"" <<
'\n';
1007 out <<
"end" <<
'\n';
1018 template <
int dim,
int spacedim>
1021 std::ostream & out)
const
1049 out <<
"@f$NOD" <<
'\n' << n_vertices <<
'\n';
1054 for (
unsigned int i = 0; i <
vertices.size(); ++i)
1059 for (
unsigned int d = spacedim + 1;
d <= 3; ++
d)
1065 out <<
"@f$ENDNOD" <<
'\n'
1076 out << cell->active_cell_index() + 1 <<
' '
1077 << cell->reference_cell().gmsh_element_type() <<
' '
1078 << cell->material_id() <<
' ' << cell->subdomain_id() <<
' '
1079 << cell->n_vertices() <<
' ';
1083 for (
const unsigned int vertex : cell->vertex_indices())
1085 if (cell->reference_cell() == ReferenceCells::get_hypercube<dim>())
1089 else if (cell->reference_cell() == ReferenceCells::get_simplex<dim>())
1090 out << cell->vertex_index(vertex) + 1 <<
' ';
1108 out <<
"@f$ENDELM\n";
1118 template <
int dim,
int spacedim>
1121 std::ostream & out)
const
1139 std::time_t time1 = std::time(
nullptr);
1140 std::tm * time = std::localtime(&time1);
1142 <<
"# This file was generated by the deal.II library." <<
'\n'
1143 <<
"# Date = " << time->tm_year + 1900 <<
"/" << time->tm_mon + 1
1144 <<
"/" << time->tm_mday <<
'\n'
1145 <<
"# Time = " << time->tm_hour <<
":" << std::setw(2) << time->tm_min
1146 <<
":" << std::setw(2) << time->tm_sec <<
'\n'
1148 <<
"# For a description of the UCD format see the AVS Developer's guide."
1154 out << n_vertices <<
' '
1164 for (
unsigned int i = 0; i <
vertices.size(); ++i)
1169 for (
unsigned int d = spacedim + 1;
d <= 3; ++
d)
1178 out << cell->active_cell_index() + 1 <<
' ' << cell->material_id() <<
' ';
1235 template <
int dim,
int spacedim>
1253 const int spacedim = 2;
1259 out <<
"#FIG 3.2\nLandscape\nCenter\nInches" << std::endl
1260 <<
"A4\n100.00\nSingle"
1263 <<
"-3" << std::endl
1264 <<
"# generated by deal.II GridOut class" << std::endl
1265 <<
"# reduce first number to scale up image" << std::endl
1266 <<
"1200 2" << std::endl;
1269 unsigned int colno = 32;
1270 out <<
"0 " << colno++ <<
" #ff0000" << std::endl;
1271 out <<
"0 " << colno++ <<
" #ff8000" << std::endl;
1272 out <<
"0 " << colno++ <<
" #ffd000" << std::endl;
1273 out <<
"0 " << colno++ <<
" #ffff00" << std::endl;
1274 out <<
"0 " << colno++ <<
" #c0ff00" << std::endl;
1275 out <<
"0 " << colno++ <<
" #80ff00" << std::endl;
1276 out <<
"0 " << colno++ <<
" #00f000" << std::endl;
1277 out <<
"0 " << colno++ <<
" #00f0c0" << std::endl;
1278 out <<
"0 " << colno++ <<
" #00f0ff" << std::endl;
1279 out <<
"0 " << colno++ <<
" #00c0ff" << std::endl;
1280 out <<
"0 " << colno++ <<
" #0080ff" << std::endl;
1281 out <<
"0 " << colno++ <<
" #0040ff" << std::endl;
1282 out <<
"0 " << colno++ <<
" #0000c0" << std::endl;
1283 out <<
"0 " << colno++ <<
" #5000ff" << std::endl;
1284 out <<
"0 " << colno++ <<
" #8000ff" << std::endl;
1285 out <<
"0 " << colno++ <<
" #b000ff" << std::endl;
1286 out <<
"0 " << colno++ <<
" #ff00ff" << std::endl;
1287 out <<
"0 " << colno++ <<
" #ff80ff" << std::endl;
1289 for (
unsigned int i = 0; i < 8; ++i)
1290 out <<
"0 " << colno++ <<
" #" << std::hex << 32 * i + 31 << 32 * i + 31
1291 << 32 * i + 31 << std::dec << std::endl;
1293 for (
unsigned int i = 1; i < 16; ++i)
1294 out <<
"0 " << colno++ <<
" #00" << std::hex << 16 * i + 15 << std::dec
1295 <<
"00" << std::endl;
1297 for (
unsigned int i = 1; i < 16; ++i)
1298 out <<
"0 " << colno++ <<
" #" << std::hex << 16 * i + 15 << 16 * i + 15
1299 << std::dec <<
"00" << std::endl;
1301 for (
unsigned int i = 1; i < 16; ++i)
1302 out <<
"0 " << colno++ <<
" #" << std::hex << 16 * i + 15 << std::dec
1303 <<
"0000" << std::endl;
1305 for (
unsigned int i = 1; i < 16; ++i)
1306 out <<
"0 " << colno++ <<
" #" << std::hex << 16 * i + 15 <<
"00"
1307 << 16 * i + 15 << std::dec << std::endl;
1309 for (
unsigned int i = 1; i < 16; ++i)
1310 out <<
"0 " << colno++ <<
" #0000" << std::hex << 16 * i + 15 << std::dec
1313 for (
unsigned int i = 1; i < 16; ++i)
1314 out <<
"0 " << colno++ <<
" #00" << std::hex << 16 * i + 15 << 16 * i + 15
1315 << std::dec << std::endl;
1337 out << cell->material_id() + 32;
1340 out << cell->level() + 8;
1343 out << cell->subdomain_id() + 32;
1346 out << cell->level_subdomain_id() + 32;
1355 (900 + cell->material_id()))
1361 << nv + 1 << std::endl;
1367 for (
unsigned int k = 0; k <= nv; ++k)
1371 for (
unsigned int d = 0; d < static_cast<unsigned int>(dim); ++
d)
1375 out <<
'\t' << ((
d == 0) ? val : -val);
1380 static const unsigned int face_reorder[4] = {2, 1, 3, 0};
1382 for (
const unsigned int f : face_reorder)
1410 for (
unsigned int k = 0;
1411 k < GeometryInfo<dim>::vertices_per_face;
1415 for (
unsigned int d = 0; d < static_cast<unsigned int>(dim);
1421 out <<
'\t' << ((
d == 0) ? val : -val);
1438 #ifdef DEAL_II_GMSH_WITH_API
1439 template <
int dim,
int spacedim>
1442 const std::string & filename)
const
1445 const std::array<int, 8> dealii_to_gmsh_type = {{15, 1, 2, 3, 4, 7, 6, 5}};
1448 const std::array<std::vector<unsigned int>, 8> dealii_to_gmsh = {
1455 {{0, 1, 2, 3, 4, 5}},
1456 {{0, 1, 3, 2, 4, 5, 7, 6}}}};
1461 std::vector<double> coords(3 *
vertices.size());
1462 std::vector<std::size_t> nodes(
vertices.size());
1468 for (
unsigned int d = 0;
d < spacedim; ++
d)
1469 coords[i * 3 +
d] = p[
d];
1481 using IdPair = std::pair<types::material_id, types::manifold_id>;
1482 std::map<IdPair, int> id_pair_to_entity_tag;
1483 std::vector<IdPair> all_pairs;
1485 std::set<IdPair> set_of_pairs;
1488 set_of_pairs.insert({cell->material_id(), cell->manifold_id()});
1489 for (
const auto &f : cell->face_iterators())
1491 (f->boundary_id() != 0 &&
1493 set_of_pairs.insert({f->boundary_id(), f->manifold_id()});
1495 for (
const auto l : cell->line_indices())
1497 const auto &f = cell->line(
l);
1499 (f->boundary_id() != 0 &&
1501 set_of_pairs.insert({f->boundary_id(), f->manifold_id()});
1504 all_pairs = {set_of_pairs.begin(), set_of_pairs.end()};
1507 for (
const auto &p : set_of_pairs)
1508 id_pair_to_entity_tag[p] = entity++;
1511 const auto n_entity_tags = id_pair_to_entity_tag.size();
1514 std::vector<std::vector<std::vector<std::size_t>>> element_ids(
1515 n_entity_tags, std::vector<std::vector<std::size_t>>(8));
1516 std::vector<std::vector<std::vector<std::size_t>>> element_nodes(
1517 n_entity_tags, std::vector<std::vector<std::size_t>>(8));
1520 std::size_t element_id = 1;
1522 const auto add_element = [&](
const auto &element,
const int &entity_tag) {
1523 const auto type = element->reference_cell();
1528 for (
const auto v : element->vertex_indices())
1529 element_nodes[entity_tag - 1][type].emplace_back(
1530 element->vertex_index(dealii_to_gmsh[type][v]) + 1);
1533 element_ids[entity_tag - 1][type].emplace_back(element_id);
1541 std::set<std::pair<int, int>> dim_entity_tag;
1543 auto maybe_add_element =
1544 [&](
const auto & element,
1546 const auto struct_dim = element->structure_dimension;
1550 const bool non_default_boundary_or_material_id =
1551 (boundary_or_material_id != 0 &&
1553 const bool non_default_manifold =
1555 if (struct_dim == dim || non_default_boundary_or_material_id ||
1556 non_default_manifold)
1558 const auto entity_tag =
1559 id_pair_to_entity_tag[{boundary_or_material_id,
manifold_id}];
1560 add_element(element, entity_tag);
1561 dim_entity_tag.insert({struct_dim, entity_tag});
1568 maybe_add_element(cell, cell->material_id());
1569 for (
const auto &face : cell->face_iterators())
1570 maybe_add_element(face, face->boundary_id());
1572 for (
const auto l : cell->line_indices())
1573 maybe_add_element(cell->line(
l), cell->line(
l)->boundary_id());
1578 gmsh::option::setNumber(
"General.Verbosity", 0);
1579 gmsh::model::add(
"Grid generated in deal.II");
1580 for (
const auto &p : dim_entity_tag)
1582 gmsh::model::addDiscreteEntity(p.first, p.second);
1583 gmsh::model::mesh::addNodes(p.first, p.second, nodes, coords);
1586 for (
unsigned int entity_tag = 0; entity_tag < n_entity_tags; ++entity_tag)
1587 for (
unsigned int t = 1; t < 8; ++t)
1589 const auto all_element_ids = element_ids[entity_tag][t];
1590 const auto all_element_nodes = element_nodes[entity_tag][t];
1591 const auto gmsh_t = dealii_to_gmsh_type[t];
1592 if (all_element_ids.size() > 0)
1593 gmsh::model::mesh::addElementsByType(entity_tag + 1,
1601 for (
const auto &it : dim_entity_tag)
1603 const auto &
d = it.first;
1604 const auto &entity_tag = it.second;
1605 const auto &
boundary_id = all_pairs[entity_tag - 1].first;
1606 const auto &
manifold_id = all_pairs[entity_tag - 1].second;
1608 std::string physical_name;
1619 std::string sep = physical_name !=
"" ?
", " :
"";
1622 sep +
"ManifoldID:" +
1624 const auto physical_tag =
1625 gmsh::model::addPhysicalGroup(
d, {entity_tag}, -1);
1626 if (physical_name !=
"")
1627 gmsh::model::setPhysicalName(
d, physical_tag, physical_name);
1631 gmsh::write(filename);
1656 const float camera_focus)
1659 cross_product_3d(camera_horizontal, camera_direction);
1662 camera_focus / ((
point - camera_position) * camera_direction);
1665 camera_position + phi * (
point - camera_position);
1667 return {(projection - camera_position - camera_focus * camera_direction) *
1669 (projection - camera_position - camera_focus * camera_direction) *
1676 template <
int dim,
int spacedim>
1679 std::ostream & )
const
1682 ExcMessage(
"Mesh output in SVG format is not implemented for anything "
1683 "other than two-dimensional meshes in two-dimensional "
1684 "space. That's because three-dimensional meshes are best "
1685 "viewed in programs that allow changing the viewpoint, "
1686 "but SVG format does not allow this: It is an inherently "
1687 "2d format, and for three-dimensional meshes would "
1688 "require choosing one, fixed viewpoint."
1690 "You probably want to output your mesh in a format such "
1691 "as VTK, VTU, or gnuplot."));
1700 unsigned int min_level, max_level;
1708 Assert(height != 0 || width != 0,
1709 ExcMessage(
"You have to set at least one of width and height"));
1711 unsigned int margin_in_percent = 0;
1713 margin_in_percent = 8;
1716 unsigned int cell_label_font_size;
1729 float x_max_perspective, x_min_perspective;
1730 float y_max_perspective, y_min_perspective;
1732 float x_dimension_perspective, y_dimension_perspective;
1736 double x_min =
tria.
begin()->vertex(0)[0];
1737 double x_max = x_min;
1738 double y_min =
tria.
begin()->vertex(0)[1];
1739 double y_max = y_min;
1741 double x_dimension, y_dimension;
1743 min_level = max_level =
tria.
begin()->level();
1746 std::set<unsigned int> materials;
1749 std::set<unsigned int> levels;
1752 std::set<unsigned int> subdomains;
1755 std::set<int> level_subdomains;
1763 for (
unsigned int vertex_index = 0; vertex_index < cell->n_vertices();
1766 if (cell->vertex(vertex_index)[0] < x_min)
1767 x_min = cell->vertex(vertex_index)[0];
1768 if (cell->vertex(vertex_index)[0] > x_max)
1769 x_max = cell->vertex(vertex_index)[0];
1771 if (cell->vertex(vertex_index)[1] < y_min)
1772 y_min = cell->vertex(vertex_index)[1];
1773 if (cell->vertex(vertex_index)[1] > y_max)
1774 y_max = cell->vertex(vertex_index)[1];
1777 if (
static_cast<unsigned int>(cell->level()) < min_level)
1778 min_level = cell->level();
1779 if (
static_cast<unsigned int>(cell->level()) > max_level)
1780 max_level = cell->level();
1782 materials.insert(cell->material_id());
1783 levels.insert(cell->level());
1784 if (cell->is_active())
1785 subdomains.insert(cell->subdomain_id() + 2);
1786 level_subdomains.insert(cell->level_subdomain_id() + 2);
1789 x_dimension = x_max - x_min;
1790 y_dimension = y_max - y_min;
1793 const unsigned int n_materials = materials.size();
1796 const unsigned int n_levels = levels.size();
1799 const unsigned int n_subdomains = subdomains.size();
1802 const unsigned int n_level_subdomains = level_subdomains.size();
1816 n = n_level_subdomains;
1825 camera_position[0] = 0;
1826 camera_position[1] = 0;
1827 camera_position[2] = 2. *
std::max(x_dimension, y_dimension);
1830 camera_direction[0] = 0;
1831 camera_direction[1] = 0;
1832 camera_direction[2] = -1;
1835 camera_horizontal[0] = 1;
1836 camera_horizontal[1] = 0;
1837 camera_horizontal[2] = 0;
1839 camera_focus = .5 *
std::max(x_dimension, y_dimension);
1845 const double angle_factor = 3.14159265 / 180.;
1848 camera_position_temp[1] =
1851 camera_position_temp[2] =
1855 camera_direction_temp[1] =
1858 camera_direction_temp[2] =
1862 camera_horizontal_temp[1] =
1865 camera_horizontal_temp[2] =
1869 camera_position[1] = camera_position_temp[1];
1870 camera_position[2] = camera_position_temp[2];
1872 camera_direction[1] = camera_direction_temp[1];
1873 camera_direction[2] = camera_direction_temp[2];
1875 camera_horizontal[1] = camera_horizontal_temp[1];
1876 camera_horizontal[2] = camera_horizontal_temp[2];
1879 camera_position_temp[0] =
1882 camera_position_temp[1] =
1886 camera_direction_temp[0] =
1889 camera_direction_temp[1] =
1893 camera_horizontal_temp[0] =
1896 camera_horizontal_temp[1] =
1900 camera_position[0] = camera_position_temp[0];
1901 camera_position[1] = camera_position_temp[1];
1903 camera_direction[0] = camera_direction_temp[0];
1904 camera_direction[1] = camera_direction_temp[1];
1906 camera_horizontal[0] = camera_horizontal_temp[0];
1907 camera_horizontal[1] = camera_horizontal_temp[1];
1910 camera_position[0] = x_min + .5 * x_dimension;
1911 camera_position[1] = y_min + .5 * y_dimension;
1913 camera_position[0] += 2. *
std::max(x_dimension, y_dimension) *
1916 camera_position[1] -= 2. *
std::max(x_dimension, y_dimension) *
1927 float min_level_min_vertex_distance = 0;
1932 (
static_cast<float>(
tria.
begin()->level()) /
1933 static_cast<float>(n_levels)) *
1934 std::max(x_dimension, y_dimension);
1937 projection_decomposition = svg_project_point(
1938 point, camera_position, camera_direction, camera_horizontal, camera_focus);
1940 x_max_perspective = projection_decomposition[0];
1941 x_min_perspective = projection_decomposition[0];
1943 y_max_perspective = projection_decomposition[1];
1944 y_min_perspective = projection_decomposition[1];
1948 point[0] = cell->vertex(0)[0];
1949 point[1] = cell->vertex(0)[1];
1956 (
static_cast<float>(cell->level()) /
static_cast<float>(n_levels)) *
1957 std::max(x_dimension, y_dimension);
1960 projection_decomposition = svg_project_point(
point,
1966 if (x_max_perspective < projection_decomposition[0])
1967 x_max_perspective = projection_decomposition[0];
1968 if (x_min_perspective > projection_decomposition[0])
1969 x_min_perspective = projection_decomposition[0];
1971 if (y_max_perspective < projection_decomposition[1])
1972 y_max_perspective = projection_decomposition[1];
1973 if (y_min_perspective > projection_decomposition[1])
1974 y_min_perspective = projection_decomposition[1];
1976 point[0] = cell->vertex(1)[0];
1977 point[1] = cell->vertex(1)[1];
1979 projection_decomposition = svg_project_point(
point,
1985 if (x_max_perspective < projection_decomposition[0])
1986 x_max_perspective = projection_decomposition[0];
1987 if (x_min_perspective > projection_decomposition[0])
1988 x_min_perspective = projection_decomposition[0];
1990 if (y_max_perspective < projection_decomposition[1])
1991 y_max_perspective = projection_decomposition[1];
1992 if (y_min_perspective > projection_decomposition[1])
1993 y_min_perspective = projection_decomposition[1];
1995 point[0] = cell->vertex(2)[0];
1996 point[1] = cell->vertex(2)[1];
1998 projection_decomposition = svg_project_point(
point,
2004 if (x_max_perspective < projection_decomposition[0])
2005 x_max_perspective = projection_decomposition[0];
2006 if (x_min_perspective > projection_decomposition[0])
2007 x_min_perspective = projection_decomposition[0];
2009 if (y_max_perspective < projection_decomposition[1])
2010 y_max_perspective = projection_decomposition[1];
2011 if (y_min_perspective > projection_decomposition[1])
2012 y_min_perspective = projection_decomposition[1];
2014 if (cell->n_vertices() == 4)
2016 point[0] = cell->vertex(3)[0];
2017 point[1] = cell->vertex(3)[1];
2019 projection_decomposition = svg_project_point(
point,
2025 if (x_max_perspective < projection_decomposition[0])
2026 x_max_perspective = projection_decomposition[0];
2027 if (x_min_perspective > projection_decomposition[0])
2028 x_min_perspective = projection_decomposition[0];
2030 if (y_max_perspective < projection_decomposition[1])
2031 y_max_perspective = projection_decomposition[1];
2032 if (y_min_perspective > projection_decomposition[1])
2033 y_min_perspective = projection_decomposition[1];
2036 if (
static_cast<unsigned int>(cell->level()) == min_level)
2037 min_level_min_vertex_distance = cell->minimum_vertex_distance();
2040 x_dimension_perspective = x_max_perspective - x_min_perspective;
2041 y_dimension_perspective = y_max_perspective - y_min_perspective;
2045 width =
static_cast<unsigned int>(
2046 .5 + height * (x_dimension_perspective / y_dimension_perspective));
2047 else if (height == 0)
2048 height =
static_cast<unsigned int>(
2049 .5 + width * (y_dimension_perspective / x_dimension_perspective));
2050 unsigned int additional_width = 0;
2052 unsigned int font_size =
2053 static_cast<unsigned int>(.5 + (height / 100.) * 1.75);
2054 cell_label_font_size =
static_cast<unsigned int>(
2056 min_level_min_vertex_distance /
std::min(x_dimension, y_dimension)));
2063 additional_width =
static_cast<unsigned int>(
2068 additional_width =
static_cast<unsigned int>(
2069 .5 + height * .175);
2081 out <<
"<svg width=\"" << width + additional_width <<
"\" height=\"" << height
2082 <<
"\" xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\">" <<
'\n'
2089 <<
" <linearGradient id=\"background_gradient\" gradientUnits=\"userSpaceOnUse\" x1=\"0\" y1=\"0\" x2=\"0\" y2=\""
2090 << height <<
"\">" <<
'\n'
2091 <<
" <stop offset=\"0\" style=\"stop-color:white\"/>" <<
'\n'
2092 <<
" <stop offset=\"1\" style=\"stop-color:lightsteelblue\"/>" <<
'\n'
2093 <<
" </linearGradient>" <<
'\n';
2099 out <<
"<!-- internal style sheet -->" <<
'\n'
2100 <<
"<style type=\"text/css\"><![CDATA[" <<
'\n';
2104 out <<
" rect.background{fill:url(#background_gradient)}" <<
'\n';
2106 out <<
" rect.background{fill:white}" <<
'\n';
2108 out <<
" rect.background{fill:none}" <<
'\n';
2111 out <<
" rect{fill:none; stroke:rgb(25,25,25); stroke-width:"
2113 <<
" text{font-family:Helvetica; text-anchor:middle; fill:rgb(25,25,25)}"
2115 <<
" line{stroke:rgb(25,25,25); stroke-width:"
2117 <<
" path{fill:none; stroke:rgb(25,25,25); stroke-width:"
2119 <<
" circle{fill:white; stroke:black; stroke-width:2}" <<
'\n'
2125 unsigned int labeling_index = 0;
2126 auto materials_it = materials.begin();
2127 auto levels_it = levels.begin();
2128 auto subdomains_it = subdomains.begin();
2129 auto level_subdomains_it = level_subdomains.begin();
2131 for (
unsigned int index = 0; index < n; ++index)
2136 h = .6 - (index / (n - 1.)) * .6;
2144 unsigned int i =
static_cast<unsigned int>(h * 6);
2146 double f = h * 6 - i;
2153 r = 255, g =
static_cast<unsigned int>(.5 + 255 * t);
2156 r =
static_cast<unsigned int>(.5 + 255 * q), g = 255;
2159 g = 255,
b =
static_cast<unsigned int>(.5 + 255 * t);
2162 g =
static_cast<unsigned int>(.5 + 255 * q),
b = 255;
2165 r =
static_cast<unsigned int>(.5 + 255 * t),
b = 255;
2168 r = 255,
b =
static_cast<unsigned int>(.5 + 255 * q);
2177 labeling_index = *materials_it++;
2180 labeling_index = *levels_it++;
2183 labeling_index = *subdomains_it++;
2186 labeling_index = *level_subdomains_it++;
2192 out <<
" path.p" << labeling_index <<
"{fill:rgb(" << r <<
',' << g
2193 <<
',' <<
b <<
"); "
2194 <<
"stroke:rgb(25,25,25); stroke-width:"
2197 out <<
" path.ps" << labeling_index <<
"{fill:rgb("
2198 <<
static_cast<unsigned int>(.5 + .75 * r) <<
','
2199 <<
static_cast<unsigned int>(.5 + .75 * g) <<
','
2200 <<
static_cast<unsigned int>(.5 + .75 *
b) <<
"); "
2201 <<
"stroke:rgb(20,20,20); stroke-width:"
2204 out <<
" rect.r" << labeling_index <<
"{fill:rgb(" << r <<
',' << g
2205 <<
',' <<
b <<
"); "
2206 <<
"stroke:rgb(25,25,25); stroke-width:"
2213 out <<
"]]></style>" <<
'\n' <<
'\n';
2216 out <<
" <rect class=\"background\" width=\"" << width <<
"\" height=\""
2217 << height <<
"\"/>" <<
'\n';
2221 unsigned int x_offset = 0;
2224 x_offset =
static_cast<unsigned int>(.5 + (height / 100.) *
2225 (margin_in_percent / 2.));
2227 x_offset =
static_cast<unsigned int>(.5 + height * .025);
2230 <<
" <text x=\"" << x_offset <<
"\" y=\""
2231 <<
static_cast<unsigned int>(.5 + height * .0525) <<
'\"'
2232 <<
" style=\"font-weight:100; fill:lightsteelblue; text-anchor:start; font-family:Courier; font-size:"
2233 <<
static_cast<unsigned int>(.5 + height * .045) <<
"px\">"
2235 <<
"</text>" <<
'\n';
2254 out <<
" <!-- cells -->" <<
'\n';
2256 for (
unsigned int level_index = min_level; level_index <= max_level;
2269 out <<
" class=\"p";
2271 if (!cell->is_active() &&
2278 out << cell->material_id();
2281 out << static_cast<unsigned int>(cell->level());
2284 if (cell->is_active())
2285 out << cell->subdomain_id() + 2;
2290 out << cell->level_subdomain_id() + 2;
2301 point[0] = cell->vertex(0)[0];
2302 point[1] = cell->vertex(0)[1];
2308 (
static_cast<float>(cell->level()) /
2309 static_cast<float>(n_levels)) *
2310 std::max(x_dimension, y_dimension);
2313 projection_decomposition = svg_project_point(
point,
2319 out << static_cast<unsigned int>(
2321 ((projection_decomposition[0] - x_min_perspective) /
2322 x_dimension_perspective) *
2323 (width - (width / 100.) * 2. * margin_in_percent) +
2324 ((width / 100.) * margin_in_percent))
2326 <<
static_cast<unsigned int>(
2327 .5 + height - (height / 100.) * margin_in_percent -
2328 ((projection_decomposition[1] - y_min_perspective) /
2329 y_dimension_perspective) *
2330 (height - (height / 100.) * 2. * margin_in_percent));
2334 point[0] = cell->vertex(1)[0];
2335 point[1] = cell->vertex(1)[1];
2337 projection_decomposition = svg_project_point(
point,
2343 out << static_cast<unsigned int>(
2345 ((projection_decomposition[0] - x_min_perspective) /
2346 x_dimension_perspective) *
2347 (width - (width / 100.) * 2. * margin_in_percent) +
2348 ((width / 100.) * margin_in_percent))
2350 <<
static_cast<unsigned int>(
2351 .5 + height - (height / 100.) * margin_in_percent -
2352 ((projection_decomposition[1] - y_min_perspective) /
2353 y_dimension_perspective) *
2354 (height - (height / 100.) * 2. * margin_in_percent));
2358 if (cell->n_vertices() == 4)
2360 point[0] = cell->vertex(3)[0];
2361 point[1] = cell->vertex(3)[1];
2363 projection_decomposition = svg_project_point(
point,
2369 out << static_cast<unsigned int>(
2371 ((projection_decomposition[0] - x_min_perspective) /
2372 x_dimension_perspective) *
2373 (width - (width / 100.) * 2. * margin_in_percent) +
2374 ((width / 100.) * margin_in_percent))
2376 <<
static_cast<unsigned int>(
2377 .5 + height - (height / 100.) * margin_in_percent -
2378 ((projection_decomposition[1] - y_min_perspective) /
2379 y_dimension_perspective) *
2380 (height - (height / 100.) * 2. * margin_in_percent));
2385 point[0] = cell->vertex(2)[0];
2386 point[1] = cell->vertex(2)[1];
2388 projection_decomposition = svg_project_point(
point,
2394 out << static_cast<unsigned int>(
2396 ((projection_decomposition[0] - x_min_perspective) /
2397 x_dimension_perspective) *
2398 (width - (width / 100.) * 2. * margin_in_percent) +
2399 ((width / 100.) * margin_in_percent))
2401 <<
static_cast<unsigned int>(
2402 .5 + height - (height / 100.) * margin_in_percent -
2403 ((projection_decomposition[1] - y_min_perspective) /
2404 y_dimension_perspective) *
2405 (height - (height / 100.) * 2. * margin_in_percent));
2409 point[0] = cell->vertex(0)[0];
2410 point[1] = cell->vertex(0)[1];
2412 projection_decomposition = svg_project_point(
point,
2418 out << static_cast<unsigned int>(
2420 ((projection_decomposition[0] - x_min_perspective) /
2421 x_dimension_perspective) *
2422 (width - (width / 100.) * 2. * margin_in_percent) +
2423 ((width / 100.) * margin_in_percent))
2425 <<
static_cast<unsigned int>(
2426 .5 + height - (height / 100.) * margin_in_percent -
2427 ((projection_decomposition[1] - y_min_perspective) /
2428 y_dimension_perspective) *
2429 (height - (height / 100.) * 2. * margin_in_percent));
2431 out <<
"\"/>" <<
'\n';
2438 point[0] = cell->center()[0];
2439 point[1] = cell->center()[1];
2445 (
static_cast<float>(cell->level()) /
2446 static_cast<float>(n_levels)) *
2447 std::max(x_dimension, y_dimension);
2450 const double distance_to_camera =
2451 std::sqrt(std::pow(
point[0] - camera_position[0], 2.) +
2452 std::pow(
point[1] - camera_position[1], 2.) +
2453 std::pow(
point[2] - camera_position[2], 2.));
2454 const double distance_factor =
2455 distance_to_camera / (2. *
std::max(x_dimension, y_dimension));
2457 projection_decomposition = svg_project_point(
point,
2463 const unsigned int font_size_this_cell =
2464 static_cast<unsigned int>(
2466 cell_label_font_size *
2467 std::pow(.5, cell->level() - 4. + 3.5 * distance_factor));
2471 <<
static_cast<unsigned int>(
2473 ((projection_decomposition[0] - x_min_perspective) /
2474 x_dimension_perspective) *
2475 (width - (width / 100.) * 2. * margin_in_percent) +
2476 ((width / 100.) * margin_in_percent))
2478 <<
static_cast<unsigned int>(
2479 .5 + height - (height / 100.) * margin_in_percent -
2480 ((projection_decomposition[1] - y_min_perspective) /
2481 y_dimension_perspective) *
2482 (height - (height / 100.) * 2. * margin_in_percent) +
2483 0.5 * font_size_this_cell)
2484 <<
"\" style=\"font-size:" << font_size_this_cell <<
"px\">";
2488 out << cell->level();
2495 out << cell->index();
2504 <<
static_cast<std::make_signed<types::material_id>::type
>(
2505 cell->material_id());
2513 if (cell->is_active())
2515 std::make_signed<types::subdomain_id>::type
>(
2516 cell->subdomain_id());
2529 <<
static_cast<std::make_signed<types::subdomain_id>::type
>(
2530 cell->level_subdomain_id());
2533 out <<
"</text>" <<
'\n';
2540 for (
auto faceIndex : cell->face_indices())
2542 if (cell->at_boundary(faceIndex))
2544 point[0] = cell->face(faceIndex)->vertex(0)[0];
2545 point[1] = cell->face(faceIndex)->vertex(0)[1];
2551 (
static_cast<float>(cell->level()) /
2552 static_cast<float>(n_levels)) *
2553 std::max(x_dimension, y_dimension);
2556 projection_decomposition =
2557 svg_project_point(
point,
2563 out <<
" <line x1=\""
2564 <<
static_cast<unsigned int>(
2566 ((projection_decomposition[0] -
2567 x_min_perspective) /
2568 x_dimension_perspective) *
2570 (width / 100.) * 2. * margin_in_percent) +
2571 ((width / 100.) * margin_in_percent))
2573 <<
static_cast<unsigned int>(
2575 (height / 100.) * margin_in_percent -
2576 ((projection_decomposition[1] -
2577 y_min_perspective) /
2578 y_dimension_perspective) *
2580 (height / 100.) * 2. * margin_in_percent));
2582 point[0] = cell->face(faceIndex)->vertex(1)[0];
2583 point[1] = cell->face(faceIndex)->vertex(1)[1];
2589 (
static_cast<float>(cell->level()) /
2590 static_cast<float>(n_levels)) *
2591 std::max(x_dimension, y_dimension);
2594 projection_decomposition =
2595 svg_project_point(
point,
2602 <<
static_cast<unsigned int>(
2604 ((projection_decomposition[0] -
2605 x_min_perspective) /
2606 x_dimension_perspective) *
2608 (width / 100.) * 2. * margin_in_percent) +
2609 ((width / 100.) * margin_in_percent))
2611 <<
static_cast<unsigned int>(
2613 (height / 100.) * margin_in_percent -
2614 ((projection_decomposition[1] -
2615 y_min_perspective) /
2616 y_dimension_perspective) *
2618 (height / 100.) * 2. * margin_in_percent))
2624 const double distance_to_camera = std::sqrt(
2625 std::pow(
point[0] - camera_position[0], 2.) +
2626 std::pow(
point[1] - camera_position[1], 2.) +
2627 std::pow(
point[2] - camera_position[2], 2.));
2628 const double distance_factor =
2629 distance_to_camera /
2630 (2. *
std::max(x_dimension, y_dimension));
2632 const unsigned int font_size_this_edge =
2633 static_cast<unsigned int>(
2634 .5 + .5 * cell_label_font_size *
2636 cell->level() - 4. +
2637 3.5 * distance_factor));
2639 point[0] = cell->face(faceIndex)->center()[0];
2640 point[1] = cell->face(faceIndex)->center()[1];
2646 (
static_cast<float>(cell->level()) /
2647 static_cast<float>(n_levels)) *
2648 std::max(x_dimension, y_dimension);
2651 projection_decomposition =
2652 svg_project_point(
point,
2658 const unsigned int xc =
static_cast<unsigned int>(
2660 ((projection_decomposition[0] - x_min_perspective) /
2661 x_dimension_perspective) *
2663 (width / 100.) * 2. * margin_in_percent) +
2664 ((width / 100.) * margin_in_percent));
2665 const unsigned int yc =
static_cast<unsigned int>(
2666 .5 + height - (height / 100.) * margin_in_percent -
2667 ((projection_decomposition[1] - y_min_perspective) /
2668 y_dimension_perspective) *
2670 (height / 100.) * 2. * margin_in_percent));
2672 out <<
" <circle cx=\"" << xc <<
"\" cy=\"" << yc
2673 <<
"\" r=\"" << font_size_this_edge <<
"\" />"
2676 out <<
" <text x=\"" << xc <<
"\" y=\"" << yc
2677 <<
"\" style=\"font-size:" << font_size_this_edge
2678 <<
"px\" dominant-baseline=\"middle\">"
2679 <<
static_cast<int>(
2680 cell->face(faceIndex)->boundary_id())
2681 <<
"</text>" <<
'\n';
2693 out <<
'\n' <<
" <!-- legend -->" <<
'\n';
2695 additional_width = 0;
2697 additional_width =
static_cast<unsigned int>(.5 + (height / 100.) * 2.5);
2705 unsigned int line_offset = 0;
2706 out <<
" <rect x=\"" << width + additional_width <<
"\" y=\""
2707 <<
static_cast<unsigned int>(.5 + (height / 100.) * margin_in_percent)
2709 <<
static_cast<unsigned int>(.5 + (height / 100.) *
2710 (40. - margin_in_percent))
2711 <<
"\" height=\"" <<
static_cast<unsigned int>(.5 + height * .215)
2714 out <<
" <text x=\""
2715 << width + additional_width +
2716 static_cast<unsigned int>(.5 + (height / 100.) * 1.25)
2718 <<
static_cast<unsigned int>(.5 +
2719 (height / 100.) * margin_in_percent +
2720 (++line_offset) * 1.5 * font_size)
2721 <<
"\" style=\"text-anchor:start; font-weight:bold; font-size:"
2722 << font_size <<
"px\">"
2724 <<
"</text>" <<
'\n';
2728 out <<
" <text x=\""
2729 << width + additional_width +
2730 static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2732 <<
static_cast<unsigned int>(.5 +
2733 (height / 100.) * margin_in_percent +
2734 (++line_offset) * 1.5 * font_size)
2735 <<
"\" style=\"text-anchor:start; font-style:oblique; font-size:"
2736 << font_size <<
"px\">"
2744 out <<
"</text>" <<
'\n';
2749 out <<
" <text x=\""
2750 << width + additional_width +
2751 static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2753 <<
static_cast<unsigned int>(.5 +
2754 (height / 100.) * margin_in_percent +
2755 (++line_offset) * 1.5 * font_size)
2756 <<
"\" style=\"text-anchor:start; font-style:oblique; font-size:"
2757 << font_size <<
"px\">"
2764 out <<
"</text>" <<
'\n';
2769 out <<
" <text x=\""
2770 << width + additional_width +
2771 static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2773 <<
static_cast<unsigned int>(.5 +
2774 (height / 100.) * margin_in_percent +
2775 (++line_offset) * 1.5 * font_size)
2776 <<
"\" style=\"text-anchor:start; font-style:oblique; font-size:"
2777 << font_size <<
"px\">"
2784 out <<
"</text>" <<
'\n';
2789 out <<
" <text x= \""
2790 << width + additional_width +
2791 static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2793 <<
static_cast<unsigned int>(.5 +
2794 (height / 100.) * margin_in_percent +
2795 (++line_offset) * 1.5 * font_size)
2796 <<
"\" style=\"text-anchor:start; font-style:oblique; font-size:"
2797 << font_size <<
"px\">"
2803 out <<
"</text>" <<
'\n';
2808 out <<
" <text x= \""
2809 << width + additional_width +
2810 static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2812 <<
static_cast<unsigned int>(.5 +
2813 (height / 100.) * margin_in_percent +
2814 (++line_offset) * 1.5 * font_size)
2815 <<
"\" style=\"text-anchor:start; font-style:oblique; font-size:"
2816 << font_size <<
"px\">"
2817 <<
"level_subdomain_id"
2818 <<
"</text>" <<
'\n';
2823 out <<
" <text x=\""
2824 << width + additional_width +
2825 static_cast<unsigned int>(.5 + (height / 100.) * 1.25)
2827 <<
static_cast<unsigned int>(.5 +
2828 (height / 100.) * margin_in_percent +
2829 (++line_offset) * 1.5 * font_size)
2830 <<
"\" style=\"text-anchor:start; font-weight:bold; font-size:"
2831 << font_size <<
"px\">"
2833 <<
"</text>" <<
'\n';
2835 out <<
" <text x= \""
2836 << width + additional_width +
2837 static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2839 <<
static_cast<unsigned int>(.5 +
2840 (height / 100.) * margin_in_percent +
2841 (++line_offset) * 1.5 * font_size)
2842 <<
"\" style=\"text-anchor:start; font-style:oblique; font-size:"
2843 << font_size <<
"px\">"
2845 <<
"</text>" <<
'\n';
2853 out <<
" <text x=\"" << width + additional_width <<
"\" y=\""
2854 <<
static_cast<unsigned int>(
2855 .5 + (height / 100.) * margin_in_percent + 13.75 * font_size)
2856 <<
"\" style=\"text-anchor:start; font-size:" << font_size <<
"px\">"
2865 out <<
'\n' <<
" <!-- colorbar -->" <<
'\n';
2867 out <<
" <text x=\"" << width + additional_width <<
"\" y=\""
2868 <<
static_cast<unsigned int>(
2869 .5 + (height / 100.) * (margin_in_percent + 29.) -
2871 <<
"\" style=\"text-anchor:start; font-weight:bold; font-size:"
2872 << font_size <<
"px\">";
2877 out <<
"material_id";
2880 out <<
"level_number";
2883 out <<
"subdomain_id";
2886 out <<
"level_subdomain_id";
2892 out <<
"</text>" <<
'\n';
2894 unsigned int element_height =
static_cast<unsigned int>(
2895 ((height / 100.) * (71. - 2. * margin_in_percent)) / n);
2896 unsigned int element_width =
2897 static_cast<unsigned int>(.5 + (height / 100.) * 2.5);
2899 int labeling_index = 0;
2900 auto materials_it = materials.begin();
2901 auto levels_it = levels.begin();
2902 auto subdomains_it = subdomains.begin();
2903 auto level_subdomains_it = level_subdomains.begin();
2905 for (
unsigned int index = 0; index < n; ++index)
2910 labeling_index = *materials_it++;
2913 labeling_index = *levels_it++;
2916 labeling_index = *subdomains_it++;
2919 labeling_index = *level_subdomains_it++;
2925 out <<
" <rect class=\"r" << labeling_index <<
"\" x=\""
2926 << width + additional_width <<
"\" y=\""
2927 <<
static_cast<unsigned int>(.5 + (height / 100.) *
2928 (margin_in_percent + 29)) +
2929 (n - index - 1) * element_height
2930 <<
"\" width=\"" << element_width <<
"\" height=\""
2931 << element_height <<
"\"/>" <<
'\n';
2933 out <<
" <text x=\""
2934 << width + additional_width + 1.5 * element_width <<
"\" y=\""
2935 <<
static_cast<unsigned int>(.5 + (height / 100.) *
2936 (margin_in_percent + 29)) +
2937 (n - index - 1 + .5) * element_height +
2938 static_cast<unsigned int>(.5 + font_size * .35)
2940 <<
" style=\"text-anchor:start; font-size:"
2941 <<
static_cast<unsigned int>(.5 + font_size) <<
"px";
2943 if (index == 0 || index == n - 1)
2944 out <<
"; font-weight:bold";
2946 out <<
"\">" << labeling_index;
2953 out <<
"</text>" <<
'\n';
2961 out <<
'\n' <<
"</svg>";
2977 template <
int dim,
int spacedim>
2980 std::ostream & out)
const
2987 const std::time_t time1 = std::time(
nullptr);
2988 const std::tm * time = std::localtime(&time1);
2992 <<
"\n# This file was generated by the deal.II library."
2993 <<
"\n# Date = " << time->tm_year + 1900 <<
"/" << std::setfill(
'0')
2994 << std::setw(2) << time->tm_mon + 1 <<
"/" << std::setfill(
'0')
2995 << std::setw(2) << time->tm_mday <<
"\n# Time = " << std::setfill(
'0')
2996 << std::setw(2) << time->tm_hour <<
":" << std::setfill(
'0')
2997 << std::setw(2) << time->tm_min <<
":" << std::setfill(
'0')
2998 << std::setw(2) << time->tm_sec <<
"\n#"
2999 <<
"\n# For a description of the MathGL script format see the MathGL manual. "
3001 <<
"\n# Note: This file is understood by MathGL v2.1 and higher only, and can "
3002 <<
"\n# be quickly viewed in a graphical environment using \'mglview\'. "
3008 const std::string axes =
"xyz";
3023 out <<
"\nsetsize 800 800";
3024 out <<
"\nrotate 0 0";
3027 out <<
"\nsetsize 800 800";
3028 out <<
"\nrotate 60 40";
3037 <<
"\n# Vertex ordering."
3038 <<
"\n# list <vertex order> <vertex indices>"
3047 out <<
"\nlist f 0 1 2 3" <<
'\n';
3051 <<
"\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"
3060 <<
"\n# List of vertices."
3061 <<
"\n# list <id> <vertices>"
3069 for (
unsigned int i = 0; i < dim; ++i)
3076 out <<
"\nlist " << axes[i] << cell->active_cell_index() <<
" ";
3078 out << cell->vertex(j)[i] <<
" ";
3085 <<
"\n# List of cells to quadplot."
3086 <<
"\n# quadplot <vertex order> <id> <style>"
3090 out <<
"\nquadplot f ";
3091 for (
unsigned int j = 0; j < dim; ++j)
3092 out << axes[j] << i <<
" ";
3117 template <
int dim,
int spacedim,
typename ITERATOR,
typename END>
3119 generate_triangulation_patches(
3125 for (; cell !=
end; ++cell)
3130 patch.
data.reinit(5, cell->n_vertices());
3132 for (
const unsigned int v : cell->vertex_indices())
3134 patch.
vertices[v] = cell->vertex(v);
3135 patch.
data(0, v) = cell->level();
3137 static_cast<std::make_signed<types::manifold_id>::type
>(
3138 cell->manifold_id());
3140 static_cast<std::make_signed<types::material_id>::type
>(
3141 cell->material_id());
3142 if (cell->is_active())
3144 static_cast<std::make_signed<types::subdomain_id>::type
>(
3145 cell->subdomain_id());
3147 patch.
data(3, v) = -1;
3149 static_cast<std::make_signed<types::subdomain_id>::type
>(
3150 cell->level_subdomain_id());
3152 patches.push_back(patch);
3158 std::vector<std::string>
3159 triangulation_patch_data_names()
3161 std::vector<std::string> v(5);
3166 v[4] =
"level_subdomain";
3173 std::vector<typename Triangulation<3, 3>::active_line_iterator>
3176 std::vector<typename Triangulation<3, 3>::active_line_iterator> res;
3178 std::vector<bool> flags;
3183 for (
const auto l : face->line_indices())
3185 const auto line = face->line(
l);
3186 if (line->user_flag_set() || line->has_children())
3189 line->set_user_flag();
3190 if (line->at_boundary())
3191 res.emplace_back(line);
3202 template <
int dim,
int spacedim>
3203 std::vector<typename Triangulation<dim, spacedim>::active_line_iterator>
3215 std::vector<typename Triangulation<3, 3>::active_line_iterator>
3218 std::vector<typename Triangulation<3, 3>::active_line_iterator> res;
3220 std::vector<bool> flags;
3225 for (
const auto l : face->line_indices())
3227 const auto line = face->line(
l);
3228 if (line->user_flag_set() || line->has_children())
3231 line->set_user_flag();
3233 (line->boundary_id() != 0 &&
3235 res.emplace_back(line);
3245 template <
int dim,
int spacedim>
3246 std::vector<typename Triangulation<dim, spacedim>::active_line_iterator>
3257 template <
int dim,
int spacedim>
3258 std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3261 std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3268 res.push_back(face);
3279 template <
int dim,
int spacedim>
3280 std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3283 std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3290 (face->boundary_id() != 0 &&
3292 res.push_back(face);
3300 template <
int dim,
int spacedim>
3303 std::ostream & out)
const
3310 const auto n_vertices =
vertices.size();
3312 out <<
"# vtk DataFile Version 3.0\n"
3313 <<
"Triangulation generated with deal.II\n"
3315 <<
"DATASET UNSTRUCTURED_GRID\n"
3316 <<
"POINTS " << n_vertices <<
" double\n";
3322 for (
unsigned int d = spacedim + 1;
d <= 3; ++
d)
3328 get_relevant_face_iterators(
tria) :
3329 get_boundary_face_iterators(
tria);
3331 get_relevant_edge_iterators(
tria) :
3332 get_boundary_edge_iterators(
tria);
3338 "At least one of the flags (output_cells, output_faces, output_edges) has to be enabled!"));
3355 cells_size += cell->n_vertices() + 1;
3358 for (
const auto &face : faces)
3359 cells_size += face->n_vertices() + 1;
3362 for (
const auto &edge : edges)
3363 cells_size += edge->n_vertices() + 1;
3367 out <<
"\nCELLS " <<
n_cells <<
' ' << cells_size <<
'\n';
3382 static const std::array<int, 8> deal_to_vtk_cell_type = {
3383 {1, 3, 5, 9, 10, 14, 13, 12}};
3384 static const std::array<unsigned int, 8> vtk_to_deal_hypercube = {
3385 {0, 1, 3, 2, 4, 5, 7, 6}};
3391 out << cell->n_vertices();
3392 for (
const unsigned int i : cell->vertex_indices())
3401 out << cell->vertex_index(vtk_to_deal_hypercube[i]);
3405 out << cell->vertex_index(i);
3408 static const std::array<unsigned int, 5> permutation_table{
3410 out << cell->vertex_index(permutation_table[i]);
3418 for (
const auto &face : faces)
3420 out << face->n_vertices();
3421 for (
const unsigned int i : face->vertex_indices())
3425 face->n_vertices() ?
3426 vtk_to_deal_hypercube[i] :
3432 for (
const auto &edge : edges)
3435 for (
const unsigned int i : edge->vertex_indices())
3436 out <<
' ' << edge->vertex_index(i);
3441 out <<
"\nCELL_TYPES " <<
n_cells <<
'\n';
3445 out << deal_to_vtk_cell_type[static_cast<int>(cell->reference_cell())]
3451 for (
const auto &face : faces)
3452 out << deal_to_vtk_cell_type[static_cast<int>(face->reference_cell())]
3458 for (
const auto &edge : edges)
3459 out << deal_to_vtk_cell_type[static_cast<int>(edge->reference_cell())]
3462 out <<
"\n\nCELL_DATA " <<
n_cells <<
'\n'
3463 <<
"SCALARS MaterialID int 1\n"
3464 <<
"LOOKUP_TABLE default\n";
3471 out << static_cast<std::make_signed<types::material_id>::type>(
3472 cell->material_id())
3479 for (
const auto &face : faces)
3481 out << static_cast<std::make_signed<types::boundary_id>::type>(
3482 face->boundary_id())
3489 for (
const auto &edge : edges)
3491 out << static_cast<std::make_signed<types::boundary_id>::type>(
3492 edge->boundary_id())
3497 out <<
"\n\nSCALARS ManifoldID int 1\n"
3498 <<
"LOOKUP_TABLE default\n";
3505 out << static_cast<std::make_signed<types::manifold_id>::type>(
3506 cell->manifold_id())
3513 for (
const auto &face : faces)
3515 out << static_cast<std::make_signed<types::manifold_id>::type>(
3516 face->manifold_id())
3523 for (
const auto &edge : edges)
3525 out << static_cast<std::make_signed<types::manifold_id>::type>(
3526 edge->manifold_id())
3539 template <
int dim,
int spacedim>
3542 std::ostream & out)
const
3550 std::vector<DataOutBase::Patch<dim, spacedim>> patches;
3557 triangulation_patch_data_names(),
3559 std::tuple<
unsigned int,
3567 out <<
" </UnstructuredGrid>\n";
3568 out <<
"<dealiiData encoding=\"base64\">";
3569 std::stringstream outstring;
3570 boost::archive::binary_oarchive ia(outstring);
3574 out <<
"\n</dealiiData>\n";
3575 out <<
"</VTKFile>\n";
3586 template <
int dim,
int spacedim>
3590 const std::string & filename_without_extension,
3591 const bool view_levels,
3592 const bool include_artificial)
const
3594 std::vector<DataOutBase::Patch<dim, spacedim>> patches;
3595 const unsigned int n_datasets = 4;
3596 std::vector<std::string> data_names;
3597 data_names.emplace_back(
"level");
3598 data_names.emplace_back(
"subdomain");
3599 data_names.emplace_back(
"level_subdomain");
3600 data_names.emplace_back(
"proc_writing");
3614 if (cell->has_children())
3616 if (!include_artificial &&
3620 else if (!include_artificial)
3622 if (cell->has_children() &&
3625 else if (cell->is_active() &&
3626 cell->level_subdomain_id() ==
3633 patch.
data.reinit(n_datasets, n_q_points);
3637 for (
unsigned int vertex = 0; vertex < n_q_points; ++vertex)
3639 patch.
vertices[vertex] = cell->vertex(vertex);
3640 patch.
data(0, vertex) = cell->level();
3641 if (cell->is_active())
3642 patch.
data(1, vertex) =
static_cast<double>(
3643 static_cast<std::make_signed<types::subdomain_id>::type
>(
3644 cell->subdomain_id()));
3646 patch.
data(1, vertex) = -1.0;
3647 patch.
data(2, vertex) =
static_cast<double>(
3648 static_cast<std::make_signed<types::subdomain_id>::type
>(
3649 cell->level_subdomain_id()));
3655 patches.push_back(patch);
3661 std::string new_file = filename_without_extension +
".vtu";
3665 new_file = filename_without_extension +
".proc" +
3670 if (tr->locally_owned_subdomain() == 0)
3672 std::vector<std::string> filenames;
3677 std::size_t pos = filename_without_extension.find_last_of(
'/');
3678 if (pos == std::string::npos)
3682 const unsigned int n_procs =
3684 for (
unsigned int i = 0; i < n_procs; ++i)
3685 filenames.push_back(filename_without_extension.substr(pos) +
3689 const std::string pvtu_filename =
3690 (filename_without_extension +
".pvtu");
3691 std::ofstream pvtu_output(pvtu_filename.c_str());
3710 std::ofstream out(new_file.c_str());
3712 std::tuple<
unsigned int,
3774 template <
int dim,
int spacedim>
3779 unsigned int n_faces = 0;
3782 if ((face->at_boundary()) && (face->boundary_id() != 0))
3790 template <
int dim,
int spacedim>
3797 std::vector<bool> line_flags;
3801 .clear_user_flags_line();
3803 unsigned int n_lines = 0;
3806 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
3807 if (cell->line(
l)->at_boundary() && (cell->line(
l)->boundary_id() != 0) &&
3808 (cell->line(
l)->user_flag_set() ==
false))
3811 cell->line(
l)->set_user_flag();
3826 const unsigned int next_element_index,
3827 std::ostream &)
const
3829 return next_element_index;
3835 const unsigned int next_element_index,
3836 std::ostream &)
const
3838 return next_element_index;
3843 const unsigned int next_element_index,
3844 std::ostream &)
const
3846 return next_element_index;
3852 const unsigned int next_element_index,
3853 std::ostream &)
const
3855 return next_element_index;
3860 const unsigned int next_element_index,
3861 std::ostream &)
const
3863 return next_element_index;
3869 const unsigned int next_element_index,
3870 std::ostream &)
const
3872 return next_element_index;
3878 const unsigned int next_element_index,
3879 std::ostream &)
const
3881 return next_element_index;
3886 const unsigned int next_element_index,
3887 std::ostream &)
const
3889 return next_element_index;
3894 template <
int dim,
int spacedim>
3897 const unsigned int next_element_index,
3898 std::ostream & out)
const
3900 unsigned int current_element_index = next_element_index;
3903 if (face->at_boundary() && (face->boundary_id() != 0))
3905 out << current_element_index <<
' '
3906 << face->reference_cell().gmsh_element_type() <<
' ';
3907 out << static_cast<unsigned int>(face->boundary_id()) <<
' '
3908 <<
static_cast<unsigned int>(face->boundary_id()) <<
' '
3909 << face->n_vertices();
3911 for (
unsigned int vertex : face->vertex_indices())
3915 << face->vertex_index(
3920 out <<
' ' << face->vertex_index(vertex) + 1;
3926 ++current_element_index;
3928 return current_element_index;
3933 template <
int dim,
int spacedim>
3936 const unsigned int next_element_index,
3937 std::ostream & out)
const
3939 unsigned int current_element_index = next_element_index;
3944 std::vector<bool> line_flags;
3948 .clear_user_flags_line();
3951 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
3952 if (cell->line(
l)->at_boundary() && (cell->line(
l)->boundary_id() != 0) &&
3953 (cell->line(
l)->user_flag_set() ==
false))
3955 out << next_element_index <<
' '
3957 out << static_cast<unsigned int>(cell->line(
l)->boundary_id()) <<
' '
3958 <<
static_cast<unsigned int>(cell->line(
l)->boundary_id())
3961 for (
unsigned int vertex = 0; vertex < 2; ++vertex)
3963 << cell->line(
l)->vertex_index(
3971 ++current_element_index;
3972 cell->line(
l)->set_user_flag();
3980 return current_element_index;
3987 const unsigned int next_element_index,
3988 std::ostream &)
const
3990 return next_element_index;
3995 const unsigned int next_element_index,
3996 std::ostream &)
const
3998 return next_element_index;
4003 const unsigned int next_element_index,
4004 std::ostream &)
const
4006 return next_element_index;
4011 const unsigned int next_element_index,
4012 std::ostream &)
const
4014 return next_element_index;
4019 const unsigned int next_element_index,
4020 std::ostream &)
const
4022 return next_element_index;
4028 const unsigned int next_element_index,
4029 std::ostream &)
const
4031 return next_element_index;
4037 const unsigned int next_element_index,
4038 std::ostream &)
const
4040 return next_element_index;
4045 const unsigned int next_element_index,
4046 std::ostream &)
const
4048 return next_element_index;
4053 template <
int dim,
int spacedim>
4056 const unsigned int next_element_index,
4057 std::ostream & out)
const
4059 unsigned int current_element_index = next_element_index;
4063 if (face->at_boundary() && (face->boundary_id() != 0))
4065 out << current_element_index <<
" "
4066 <<
static_cast<unsigned int>(face->boundary_id()) <<
" ";
4079 for (
unsigned int vertex = 0;
4080 vertex < GeometryInfo<dim>::vertices_per_face;
4082 out << face->vertex_index(
4088 ++current_element_index;
4090 return current_element_index;
4095 template <
int dim,
int spacedim>
4098 const unsigned int next_element_index,
4099 std::ostream & out)
const
4101 unsigned int current_element_index = next_element_index;
4106 std::vector<bool> line_flags;
4110 .clear_user_flags_line();
4113 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
4114 if (cell->line(
l)->at_boundary() && (cell->line(
l)->boundary_id() != 0) &&
4115 (cell->line(
l)->user_flag_set() ==
false))
4117 out << current_element_index <<
" "
4118 <<
static_cast<unsigned int>(cell->line(
l)->boundary_id())
4121 for (
unsigned int vertex = 0; vertex < 2; ++vertex)
4122 out << cell->line(
l)->vertex_index(
4131 ++current_element_index;
4132 cell->line(
l)->set_user_flag();
4139 return current_element_index;
4156 template <
int spacedim>
4160 while (points.size() > 2)
4163 first_difference /= first_difference.
norm();
4165 second_difference /= second_difference.
norm();
4167 if ((first_difference - second_difference).
norm() < 1
e-10)
4168 points.erase(points.begin() + 1);
4176 template <
int spacedim>
4188 out <<
"# cell " << cell <<
'\n';
4190 out << cell->vertex(0) <<
' ' << cell->level() <<
' '
4191 << cell->material_id() <<
'\n'
4192 << cell->vertex(1) <<
' ' << cell->level() <<
' '
4193 << cell->material_id() <<
'\n'
4206 template <
int spacedim>
4217 const unsigned int n_additional_points =
4219 const unsigned int n_points = 2 + n_additional_points;
4224 std::vector<
Point<dim - 1>> boundary_points;
4225 if (mapping !=
nullptr)
4227 boundary_points.resize(n_points);
4228 boundary_points[0][0] = 0;
4229 boundary_points[n_points - 1][0] = 1;
4230 for (
unsigned int i = 1; i < n_points - 1; ++i)
4231 boundary_points[i](0) = 1. * i / (n_points - 1);
4233 std::vector<double> dummy_weights(n_points, 1. / n_points);
4234 Quadrature<dim - 1> quadrature(boundary_points, dummy_weights);
4237 ReferenceCells::get_hypercube<dim>(), quadrature);
4243 out <<
"# cell " << cell <<
'\n';
4245 if (mapping ==
nullptr ||
4257 << cell->level() <<
' ' << cell->material_id() <<
'\n';
4258 out << cell->vertex(0) <<
' ' << cell->level() <<
' '
4259 << cell->material_id() <<
'\n'
4267 for (
const unsigned int face_no :
4270 const typename ::Triangulation<dim,
4271 spacedim>::face_iterator
4272 face = cell->face(face_no);
4273 if (dim != spacedim || face->at_boundary() ||
4279 std::vector<Point<spacedim>> line_points;
4282 const unsigned int offset = face_no * n_points;
4283 for (
unsigned int i = 0; i < n_points; ++i)
4284 line_points.push_back(
4286 cell, q_projector.
point(offset + i)));
4287 internal::remove_colinear_points(line_points);
4290 out <<
point <<
' ' << cell->level() <<
' '
4291 << cell->material_id() <<
'\n';
4293 out <<
'\n' <<
'\n';
4299 out << face->vertex(0) <<
' ' << cell->level() <<
' '
4300 << cell->material_id() <<
'\n'
4301 << face->vertex(1) <<
' ' << cell->level() <<
' '
4302 << cell->material_id() <<
'\n'
4318 template <
int spacedim>
4329 const unsigned int n_additional_points =
4331 const unsigned int n_points = 2 + n_additional_points;
4335 std::unique_ptr<Quadrature<dim>> q_projector;
4336 std::vector<Point<1>> boundary_points;
4337 if (mapping !=
nullptr)
4339 boundary_points.resize(n_points);
4340 boundary_points[0][0] = 0;
4341 boundary_points[n_points - 1][0] = 1;
4342 for (
unsigned int i = 1; i < n_points - 1; ++i)
4343 boundary_points[i](0) = 1. * i / (n_points - 1);
4345 std::vector<double> dummy_weights(n_points, 1. / n_points);
4349 QIterated<dim - 1> quadrature(quadrature1d, 1);
4350 q_projector = std::make_unique<Quadrature<dim>>(
4352 ReferenceCells::get_hypercube<dim>(), quadrature));
4358 out <<
"# cell " << cell <<
'\n';
4360 if (mapping ==
nullptr || n_points == 2 ||
4361 (!cell->has_boundary_lines() &&
4367 out << cell->vertex(0) <<
' ' << cell->level() <<
' '
4368 << cell->material_id() <<
'\n'
4369 << cell->vertex(1) <<
' ' << cell->level() <<
' '
4370 << cell->material_id() <<
'\n'
4371 << cell->vertex(5) <<
' ' << cell->level() <<
' '
4372 << cell->material_id() <<
'\n'
4373 << cell->vertex(4) <<
' ' << cell->level() <<
' '
4374 << cell->material_id() <<
'\n'
4375 << cell->vertex(0) <<
' ' << cell->level() <<
' '
4376 << cell->material_id() <<
'\n'
4379 out << cell->vertex(2) <<
' ' << cell->level() <<
' '
4380 << cell->material_id() <<
'\n'
4381 << cell->vertex(3) <<
' ' << cell->level() <<
' '
4382 << cell->material_id() <<
'\n'
4383 << cell->vertex(7) <<
' ' << cell->level() <<
' '
4384 << cell->material_id() <<
'\n'
4385 << cell->vertex(6) <<
' ' << cell->level() <<
' '
4386 << cell->material_id() <<
'\n'
4387 << cell->vertex(2) <<
' ' << cell->level() <<
' '
4388 << cell->material_id() <<
'\n'
4392 out << cell->vertex(0) <<
' ' << cell->level() <<
' '
4393 << cell->material_id() <<
'\n'
4394 << cell->vertex(2) <<
' ' << cell->level() <<
' '
4395 << cell->material_id() <<
'\n'
4397 out << cell->vertex(1) <<
' ' << cell->level() <<
' '
4398 << cell->material_id() <<
'\n'
4399 << cell->vertex(3) <<
' ' << cell->level() <<
' '
4400 << cell->material_id() <<
'\n'
4402 out << cell->vertex(5) <<
' ' << cell->level() <<
' '
4403 << cell->material_id() <<
'\n'
4404 << cell->vertex(7) <<
' ' << cell->level() <<
' '
4405 << cell->material_id() <<
'\n'
4407 out << cell->vertex(4) <<
' ' << cell->level() <<
' '
4408 << cell->material_id() <<
'\n'
4409 << cell->vertex(6) <<
' ' << cell->level() <<
' '
4410 << cell->material_id() <<
'\n'
4416 for (
const unsigned int v : {0, 1, 2, 0, 3, 2})
4418 out << cell->vertex(v) <<
' ' << cell->level() <<
' '
4419 << cell->material_id() <<
'\n';
4423 for (
const unsigned int v : {3, 1})
4425 out << cell->vertex(v) <<
' ' << cell->level() <<
' '
4426 << cell->material_id() <<
'\n';
4437 for (
const unsigned int v : {0, 1, 2, 0, 3, 4, 5, 3})
4439 out << cell->vertex(v) <<
' ' << cell->level() <<
' '
4440 << cell->material_id() <<
'\n';
4444 for (
const unsigned int v : {1, 4})
4446 out << cell->vertex(v) <<
' ' << cell->level() <<
' '
4447 << cell->material_id() <<
'\n';
4451 for (
const unsigned int v : {2, 5})
4453 out << cell->vertex(v) <<
' ' << cell->level() <<
' '
4454 << cell->material_id() <<
'\n';
4461 for (
const unsigned int v : {0, 1, 3, 2, 0, 4, 1})
4463 out << cell->vertex(v) <<
' ' << cell->level() <<
' '
4464 << cell->material_id() <<
'\n';
4468 for (
const unsigned int v : {2, 4, 3})
4470 out << cell->vertex(v) <<
' ' << cell->level() <<
' '
4471 << cell->material_id() <<
'\n';
4482 for (
const unsigned int face_no :
4485 const typename ::Triangulation<dim,
4486 spacedim>::face_iterator
4487 face = cell->face(face_no);
4489 if (face->at_boundary() &&
4492 const unsigned int offset = face_no * n_points * n_points;
4493 for (
unsigned int i = 0; i < n_points - 1; ++i)
4494 for (
unsigned int j = 0; j < n_points - 1; ++j)
4499 q_projector->point(offset + i * n_points + j));
4500 out << p0 <<
' ' << cell->level() <<
' '
4501 << cell->material_id() <<
'\n';
4505 offset + (i + 1) * n_points + j)))
4506 <<
' ' << cell->level() <<
' '
4507 << cell->material_id() <<
'\n';
4511 offset + (i + 1) * n_points + j + 1)))
4512 <<
' ' << cell->level() <<
' '
4513 << cell->material_id() <<
'\n';
4516 q_projector->point(offset + i * n_points +
4518 <<
' ' << cell->level() <<
' '
4519 << cell->material_id() <<
'\n';
4521 out << p0 <<
' ' << cell->level() <<
' '
4522 << cell->material_id() <<
'\n';
4523 out <<
'\n' <<
'\n';
4528 for (
unsigned int l = 0;
4529 l < GeometryInfo<dim>::lines_per_face;
4532 const typename ::Triangulation<dim, spacedim>::
4533 line_iterator line = face->line(
l);
4536 &
v1 = line->vertex(1);
4537 if (line->at_boundary() ||
4543 std::vector<Point<spacedim>> line_points;
4552 for (
unsigned int i = 0; i < n_points; ++i)
4553 line_points.push_back(
4556 (1 - boundary_points[i][0]) * u0 +
4557 boundary_points[i][0] * u1));
4558 internal::remove_colinear_points(line_points);
4560 out <<
point <<
' ' << cell->level() <<
' '
4561 <<
static_cast<unsigned int>(
4562 cell->material_id())
4566 out <<
v0 <<
' ' << cell->level() <<
' '
4567 << cell->material_id() <<
'\n'
4568 <<
v1 <<
' ' << cell->level() <<
' '
4569 << cell->material_id() <<
'\n';
4571 out <<
'\n' <<
'\n';
4588 template <
int dim,
int spacedim>
4612 const unsigned int l)
4632 write_eps(const ::Triangulation<1, 2> &,
4642 write_eps(const ::Triangulation<1, 3> &,
4652 write_eps(const ::Triangulation<2, 3> &,
4663 template <
int dim,
int spacedim>
4671 using LineList = std::list<LineEntry>;
4712 for (
const unsigned int line_no : cell->line_indices())
4714 typename ::Triangulation<dim, spacedim>::line_iterator
4715 line = cell->line(line_no);
4727 if (!line->has_children() &&
4728 (mapping ==
nullptr || !line->at_boundary()))
4748 line_list.emplace_back(
4749 Point<2>(line->vertex(0)(0), line->vertex(0)(1)),
4750 Point<2>(line->vertex(1)(0), line->vertex(1)(1)),
4751 line->user_flag_set(),
4761 if (mapping !=
nullptr)
4768 std::vector<
Point<dim - 1>> boundary_points(n_points);
4770 for (
unsigned int i = 0; i < n_points; ++i)
4771 boundary_points[i](0) = 1. * (i + 1) / (n_points + 1);
4773 Quadrature<dim - 1> quadrature(boundary_points);
4776 ReferenceCells::get_hypercube<dim>(), quadrature));
4783 for (
const unsigned int face_no :
4786 const typename ::Triangulation<dim, spacedim>::
4787 face_iterator face = cell->face(face_no);
4789 if (face->at_boundary())
4799 const unsigned int offset = face_no * n_points;
4800 for (
unsigned int i = 0; i < n_points; ++i)
4804 cell, q_projector.point(offset + i)));
4805 const Point<2> p1(p1_dim(0), p1_dim(1));
4807 line_list.emplace_back(p0,
4809 face->user_flag_set(),
4816 const Point<2> p1(p1_dim(0), p1_dim(1));
4817 line_list.emplace_back(p0,
4819 face->user_flag_set(),
4854 const double turn_angle = eps_flags_3.
turn_angle;
4856 -std::sin(z_angle * 2. * pi / 360.) *
4857 std::sin(turn_angle * 2. * pi / 360.),
4858 +std::sin(z_angle * 2. * pi / 360.) *
4859 std::cos(turn_angle * 2. * pi / 360.),
4860 -std::cos(z_angle * 2. * pi / 360.));
4868 ((
Point<dim>(0, 0, 1) * view_direction) * view_direction);
4877 ((
Point<dim>(1, 0, 0) * view_direction) * view_direction) -
4878 ((
Point<dim>(1, 0, 0) * unit_vector1) * unit_vector1));
4883 for (
const unsigned int line_no : cell->line_indices())
4885 typename ::Triangulation<dim, spacedim>::line_iterator
4886 line = cell->line(line_no);
4887 line_list.emplace_back(
4888 Point<2>(line->vertex(0) * unit_vector2,
4889 line->vertex(0) * unit_vector1),
4890 Point<2>(line->vertex(1) * unit_vector2,
4891 line->vertex(1) * unit_vector1),
4892 line->user_flag_set(),
4909 double x_max = x_min;
4911 double y_max = y_min;
4912 unsigned int max_level = line_list.begin()->level;
4914 for (LineList::const_iterator line = line_list.begin();
4915 line != line_list.end();
4918 x_min =
std::min(x_min, line->first(0));
4919 x_min =
std::min(x_min, line->second(0));
4921 x_max =
std::max(x_max, line->first(0));
4922 x_max =
std::max(x_max, line->second(0));
4924 y_min =
std::min(y_min, line->first(1));
4925 y_min =
std::min(y_min, line->second(1));
4927 y_max =
std::max(y_max, line->first(1));
4928 y_max =
std::max(y_max, line->second(1));
4930 max_level =
std::max(max_level, line->level);
4938 const double scale =
4939 (eps_flags_base.
size /
4950 std::time_t time1 = std::time(
nullptr);
4951 std::tm * time = std::localtime(&time1);
4952 out <<
"%!PS-Adobe-2.0 EPSF-1.2" <<
'\n'
4953 <<
"%%Title: deal.II Output" <<
'\n'
4954 <<
"%%Creator: the deal.II library" <<
'\n'
4955 <<
"%%Creation Date: " << time->tm_year + 1900 <<
"/"
4956 << time->tm_mon + 1 <<
"/" << time->tm_mday <<
" - "
4957 << time->tm_hour <<
":" << std::setw(2) << time->tm_min <<
":"
4958 << std::setw(2) << time->tm_sec <<
'\n'
4959 <<
"%%BoundingBox: "
4963 <<
static_cast<unsigned int>(
4966 <<
static_cast<unsigned int>(
4976 out <<
"/m {moveto} bind def" <<
'\n'
4977 <<
"/x {lineto stroke} bind def" <<
'\n'
4978 <<
"/b {0 0 0 setrgbcolor} def" <<
'\n'
4979 <<
"/r {1 0 0 setrgbcolor} def" <<
'\n';
4986 out <<
"/l { neg " << (max_level) <<
" add "
4987 << (0.66666 /
std::max(1U, (max_level - 1)))
4988 <<
" mul 1 0.8 sethsbcolor} def" <<
'\n';
5002 << (
"/R {rmoveto} bind def\n"
5003 "/Symbol-Oblique /Symbol findfont [1 0 .167 1 0 0] makefont\n"
5004 "dup length dict begin {1 index /FID eq {pop pop} {def} ifelse} forall\n"
5005 "currentdict end definefont\n"
5006 "/MFshow {{dup dup 0 get findfont exch 1 get scalefont setfont\n"
5007 "[ currentpoint ] exch dup 2 get 0 exch rmoveto dup dup 5 get exch 4 get\n"
5008 "{show} {stringwidth pop 0 rmoveto}ifelse dup 3 get\n"
5009 "{2 get neg 0 exch rmoveto pop} {pop aload pop moveto}ifelse} forall} bind def\n"
5010 "/MFwidth {0 exch {dup 3 get{dup dup 0 get findfont exch 1 get scalefont setfont\n"
5011 "5 get stringwidth pop add}\n"
5012 "{pop} ifelse} forall} bind def\n"
5013 "/MCshow { currentpoint stroke m\n"
5014 "exch dup MFwidth -2 div 3 -1 roll R MFshow } def\n")
5018 out <<
"%%EndProlog" <<
'\n' <<
'\n';
5021 out << eps_flags_base.
line_width <<
" setlinewidth" <<
'\n';
5025 const Point<2> offset(x_min, y_min);
5027 for (LineList::const_iterator line = line_list.begin();
5028 line != line_list.end();
5031 out << line->level <<
" l " << (line->first - offset) *
scale <<
" m "
5032 << (line->second - offset) *
scale <<
" x" <<
'\n';
5037 << (line->first - offset) *
scale <<
" m "
5038 << (line->second - offset) *
scale <<
" x" <<
'\n';
5044 out <<
"(Helvetica) findfont 140 scalefont setfont" <<
'\n';
5048 out << (cell->center()(0) - offset(0)) *
scale <<
' '
5049 << (cell->center()(1) - offset(1)) *
scale <<
" m" <<
'\n'
5050 <<
"[ [(Helvetica) 12.0 0.0 true true (";
5054 out << cell->index();
5057 <<
"] -6 MCshow" <<
'\n';
5064 out <<
"(Helvetica) findfont 140 scalefont setfont" <<
'\n';
5070 std::set<unsigned int> treated_vertices;
5072 for (
const unsigned int vertex_no : cell->vertex_indices())
5073 if (treated_vertices.find(cell->vertex_index(vertex_no)) ==
5074 treated_vertices.end())
5076 treated_vertices.insert(cell->vertex_index(vertex_no));
5078 out << (cell->vertex(vertex_no)(0) - offset(0)) *
scale <<
' '
5079 << (cell->vertex(vertex_no)(1) - offset(1)) *
scale
5081 <<
"[ [(Helvetica) 10.0 0.0 true true ("
5082 << cell->vertex_index(vertex_no) <<
")] "
5083 <<
"] -6 MCshow" <<
'\n';
5087 out <<
"showpage" <<
'\n';
5099 template <
int dim,
int spacedim>
5109 template <
int dim,
int spacedim>
5116 switch (output_format)
5166 template <
int dim,
int spacedim>
5177 #include "grid_out.inst"
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)
unsigned int n_boundary_faces(const Triangulation< dim, spacedim > &tria) const
GridOutFlags::Vtu vtu_flags
GridOutFlags::Eps< 2 > eps_flags_2
unsigned int write_ucd_lines(const Triangulation< dim, spacedim > &tria, const unsigned int next_element_index, std::ostream &out) const
void parse_parameters(ParameterHandler ¶m)
void write_svg(const Triangulation< 2, 2 > &tria, std::ostream &out) const
void set_flags(const GridOutFlags::DX &flags)
GridOutFlags::Eps< 1 > eps_flags_1
GridOutFlags::XFig xfig_flags
unsigned int write_ucd_faces(const Triangulation< dim, spacedim > &tria, const unsigned int next_element_index, std::ostream &out) const
void write_mathgl(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
std::string default_suffix() const
static void declare_parameters(ParameterHandler ¶m)
static OutputFormat parse_output_format(const std::string &format_name)
void write_vtk(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
GridOutFlags::Gnuplot gnuplot_flags
void write_msh(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
unsigned int write_msh_lines(const Triangulation< dim, spacedim > &tria, const unsigned int next_element_index, std::ostream &out) const
void write_eps(const Triangulation< dim, spacedim > &tria, std::ostream &out, const Mapping< dim, spacedim > *mapping=nullptr) const
static std::string get_output_format_names()
GridOutFlags::Eps< 3 > eps_flags_3
void write(const Triangulation< dim, spacedim > &tria, std::ostream &out, const OutputFormat output_format, const Mapping< dim, spacedim > *mapping=nullptr) const
void write_vtu(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
GridOutFlags::Ucd ucd_flags
void write_dx(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
std::size_t memory_consumption() const
unsigned int n_boundary_lines(const Triangulation< dim, spacedim > &tria) const
GridOutFlags::DX dx_flags
unsigned int write_msh_faces(const Triangulation< dim, spacedim > &tria, const unsigned int next_element_index, std::ostream &out) const
void write_ucd(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
GridOutFlags::Svg svg_flags
OutputFormat default_format
void write_xfig(const Triangulation< dim, spacedim > &tria, std::ostream &out, const Mapping< dim, spacedim > *mapping=nullptr) const
void write_gnuplot(const Triangulation< dim, spacedim > &tria, std::ostream &out, const Mapping< dim, spacedim > *mapping=nullptr) const
GridOutFlags::Vtk vtk_flags
GridOutFlags::Msh msh_flags
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
@ vtk
write() calls write_vtk()
@ eps
write() calls write_eps()
@ msh
write() calls write_msh()
@ xfig
write() calls write_xfig()
@ dx
write() calls write_dx()
@ ucd
write() calls write_ucd()
@ gnuplot
write() calls write_gnuplot()
@ mathgl
write() calls write_mathgl()
@ svg
write() calls write_svg()
@ none
Do nothing in write()
@ vtu
write() calls write_vtu()
GridOutFlags::MathGL mathgl_flags
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const =0
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &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)
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
const std::vector< bool > & get_used_vertices() const
void save_user_flags_line(std::ostream &out) const
void save(Archive &ar, const unsigned int version) const
cell_iterator begin(const unsigned int level=0) const
virtual types::subdomain_id locally_owned_subdomain() const
unsigned int n_active_cells() const
unsigned int n_used_vertices() const
cell_iterator end() const
const std::vector< ReferenceCell > & get_reference_cells() const
const std::vector< Point< spacedim > > & get_vertices() const
active_cell_iterator begin_active(const unsigned int level=0) const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
IteratorRange< active_cell_iterator > active_cell_iterators() const
IteratorRange< active_face_iterator > active_face_iterators() const
IteratorRange< cell_iterator > cell_iterators() const
IteratorRange< cell_iterator > cell_iterators_on_level(const unsigned int level) const
static ::ExceptionBase & ExcInvalidState()
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
DataComponentInterpretation
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_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 write_vtu_header(std::ostream &out, const VtkFlags &flags)
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_footer(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_cells(const std::vector< Patch< dim, spacedim >> &patches, StreamType &out)
Expression floor(const Expression &x)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
types::global_dof_index size_type
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
constexpr const ReferenceCell Tetrahedron
constexpr const ReferenceCell Quadrilateral
constexpr const ReferenceCell Wedge
constexpr const ReferenceCell Pyramid
constexpr const ReferenceCell Triangle
constexpr const ReferenceCell Hexahedron
constexpr const ReferenceCell Vertex
constexpr const ReferenceCell Line
VectorType::value_type * end(VectorType &V)
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
std::string encode_base64(const std::vector< unsigned char > &binary_input)
std::string compress(const std::string &input)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
const types::boundary_id invalid_boundary_id
static constexpr double PI
const types::boundary_id internal_face_boundary_id
const types::subdomain_id artificial_subdomain_id
static const unsigned int invalid_unsigned_int
const types::manifold_id flat_manifold_id
static void declare_parameters(ParameterHandler &prm)
void parse_parameters(const ParameterHandler &prm)
ReferenceCell reference_cell
unsigned int n_subdivisions
std::array< Point< spacedim >, GeometryInfo< dim >::vertices_per_cell > vertices
bool points_are_available
std::array< unsigned int, GeometryInfo< dim >::faces_per_cell > neighbors
static void declare_parameters(ParameterHandler ¶m)
void parse_parameters(ParameterHandler ¶m)
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)
static void declare_parameters(ParameterHandler ¶m)
bool color_lines_on_user_flag
unsigned int n_boundary_face_points
void parse_parameters(ParameterHandler ¶m)
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)
void parse_parameters(ParameterHandler ¶m)
bool write_cell_number_level
void parse_parameters(ParameterHandler ¶m)
bool write_vertex_numbers
void parse_parameters(ParameterHandler ¶m)
unsigned int n_extra_curved_line_points
void parse_parameters(ParameterHandler ¶m)
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)
static void declare_parameters(ParameterHandler ¶m)
bool write_additional_boundary_lines
void parse_parameters(ParameterHandler ¶m)
static void declare_parameters(ParameterHandler ¶m)
Msh(const bool write_faces=false, const bool write_lines=false)
void parse_parameters(ParameterHandler ¶m)
static void declare_parameters(ParameterHandler ¶m)
bool label_level_subdomain_id
unsigned int line_thickness
bool convert_level_number_to_height
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)
@ level_subdomain_id
Convert the level subdomain id into the cell color.
@ subdomain_id
Convert the subdomain id into the cell color.
@ material_id
Convert the material id into the cell color (default)
@ level_number
Convert the level number into the cell color.
unsigned int boundary_line_thickness
float level_height_factor
Ucd(const bool write_preamble=false, const bool write_faces=false, const bool write_lines=false)
static void declare_parameters(ParameterHandler ¶m)
void parse_parameters(ParameterHandler ¶m)
bool output_only_relevant
bool serialize_triangulation
unsigned int n_boundary_face_points
@ level_number
Convert the level into the cell color.
@ material_id
Convert the material id into the cell color.
@ level_subdomain_id
Convert the level subdomain id into the cell color.
@ subdomain_id
Convert the global subdomain id into the cell color.
void parse_parameters(ParameterHandler ¶m)
static void declare_parameters(ParameterHandler ¶m)
enum GridOutFlags::XFig::Coloring color_by
const ::Triangulation< dim, spacedim > & tria