44#ifdef DEAL_II_WITH_ZLIB
48#ifdef DEAL_II_WITH_HDF5
52#ifdef DEAL_II_WITH_NETCDF
54# include <netcdf_meta.h>
55# include <netcdf_par.h>
58#include <boost/archive/iterators/base64_from_binary.hpp>
59#include <boost/archive/iterators/transform_width.hpp>
60#include <boost/iostreams/copy.hpp>
61#include <boost/iostreams/device/back_inserter.hpp>
62#include <boost/iostreams/filtering_stream.hpp>
63#ifdef DEAL_II_WITH_ZLIB
64# include <boost/iostreams/filter/zlib.hpp>
79 <<
"Unexpected input: expected line\n <" << arg1
80 <<
">\nbut got\n <" << arg2 <<
">");
82# ifdef DEAL_II_WITH_ZLIB
83 constexpr bool deal_ii_with_zlib =
true;
85 constexpr bool deal_ii_with_zlib =
false;
89# ifdef DEAL_II_WITH_ZLIB
100 return Z_NO_COMPRESSION;
104 return Z_BEST_COMPRESSION;
106 return Z_DEFAULT_COMPRESSION;
109 return Z_NO_COMPRESSION;
113# ifdef DEAL_II_WITH_MPI
124 return boost::iostreams::zlib::no_compression;
126 return boost::iostreams::zlib::best_speed;
128 return boost::iostreams::zlib::best_compression;
130 return boost::iostreams::zlib::default_compression;
133 return boost::iostreams::zlib::no_compression;
143 template <
typename T>
145 compress_array(
const std::vector<T> &
data,
147 std::ostringstream &output)
149# ifdef DEAL_II_WITH_ZLIB
150 if (
data.size() != 0)
152 const std::size_t uncompressed_size = (
data.size() *
sizeof(
T));
160 std::numeric_limits<std::uint32_t>::max(),
164 auto compressed_data_length = compressBound(uncompressed_size);
166 std::numeric_limits<std::uint32_t>::max(),
169 std::vector<unsigned char> compressed_data(compressed_data_length);
171 int err = compress2(&compressed_data[0],
172 &compressed_data_length,
173 reinterpret_cast<const Bytef *
>(
data.data()),
175 get_zlib_compression_level(compression_level));
179 compressed_data.resize(compressed_data_length);
182 const std::uint32_t compression_header[4] = {
184 static_cast<std::uint32_t
>(uncompressed_size),
185 static_cast<std::uint32_t
>(
187 static_cast<std::uint32_t
>(
188 compressed_data_length)};
191 const std::array<std::string, 3> paddings{{
"",
"==",
"="}};
195 static_assert(
sizeof(compressed_data[0]) == 1);
196 using namespace boost::archive::iterators;
198 base64_from_binary<transform_width<const unsigned char *, 6, 8>>;
200 auto char_begin =
reinterpret_cast<const unsigned char *
>(
201 std::begin(compression_header)),
202 char_end =
reinterpret_cast<const unsigned char *
>(
203 std::end(compression_header));
205 for (
auto it =
begin; it !=
end; ++it)
208 output << paddings[(char_end - char_begin) % 3];
215 auto it =
iterator(compressed_data.data()),
216 end =
iterator(compressed_data.data() + compressed_data.size());
217 std::array<char, 128> buffer;
220 std::streamsize count = 0;
221 for (; count < static_cast<std::streamsize>(buffer.size());
229 output.write(buffer.data(), count);
232 output << paddings[compressed_data.size() % 3];
237 (void)compression_level;
240 ExcMessage(
"This function can only be called if cmake found "
241 "a working libz installation."));
255 template <
typename T>
257 vtu_stringize_array(
const std::vector<T> &
data,
260 std::ostringstream &output)
262 if (deal_ii_with_zlib &&
266 compress_array(
data, compression_level, output);
270 const auto old_precision = output.precision(precision);
271 for (
const T &el :
data)
273 output.precision(old_precision);
286 struct ParallelIntermediateHeader
289 std::uint64_t version;
290 std::uint64_t compression;
291 std::uint64_t dimension;
292 std::uint64_t space_dimension;
293 std::uint64_t n_ranks;
294 std::uint64_t n_patches;
344 SvgCell::operator<(
const SvgCell &e)
const
347 return depth >
e.depth;
385 EpsCell2d::operator<(
const EpsCell2d &e)
const
388 return depth > e.depth;
404 template <
int dim,
int spacedim,
typename Number =
double>
405 std::unique_ptr<Table<2, Number>>
406 create_global_data_table(
const std::vector<Patch<dim, spacedim>> &patches)
410 return std::make_unique<Table<2, Number>>();
419 const unsigned int n_data_sets = patches[0].points_are_available ?
420 (patches[0].data.n_rows() - spacedim) :
421 patches[0].
data.n_rows();
422 const unsigned int n_data_points =
423 std::accumulate(patches.begin(),
426 [](
const unsigned int count,
427 const Patch<dim, spacedim> &patch) {
428 return count + patch.data.n_cols();
431 std::unique_ptr<Table<2, Number>> global_data_table =
432 std::make_unique<Table<2, Number>>(n_data_sets, n_data_points);
435 unsigned int next_value = 0;
436 for (
const auto &patch : patches)
438 const unsigned int n_subdivisions = patch.n_subdivisions;
439 (void)n_subdivisions;
441 Assert((patch.data.n_rows() == n_data_sets &&
442 !patch.points_are_available) ||
443 (patch.data.n_rows() == n_data_sets + spacedim &&
444 patch.points_are_available),
446 (n_data_sets + spacedim) :
448 patch.
data.n_rows()));
449 Assert(patch.reference_cell != ReferenceCells::get_hypercube<dim>() ||
450 (n_data_sets == 0) ||
451 (patch.data.n_cols() ==
452 Utilities::fixed_power<dim>(n_subdivisions + 1)),
454 n_subdivisions + 1));
456 for (
unsigned int i = 0; i < patch.data.n_cols(); ++i, ++next_value)
457 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
458 (*global_data_table)[data_set][next_value] =
459 patch.data(data_set, i);
463 return global_data_table;
473 , node_dim(
numbers::invalid_unsigned_int)
481 , node_dim(
numbers::invalid_unsigned_int)
494 for (
unsigned int d = 0; d < dim; ++d)
498 unsigned int internal_ind;
509 internal_ind = it->second;
519 const unsigned int pt_index)
537 for (
unsigned int d = 0; d <
node_dim; ++d)
538 node_data[
node_dim * existing_point.second + d] =
539 existing_point.first[d];
547 std::vector<unsigned int> &cell_data)
const
553 cell_data[filtered_cell.first] =
554 filtered_cell.second + local_node_offset;
623 const unsigned int start,
624 const std::array<unsigned int, dim> &offsets)
628 const unsigned int base_entry =
641 const unsigned int d1 = offsets[0];
650 const unsigned int d1 = offsets[0];
651 const unsigned int d2 = offsets[1];
662 const unsigned int d1 = offsets[0];
663 const unsigned int d2 = offsets[1];
664 const unsigned int d3 = offsets[2];
687 const unsigned int start,
688 const unsigned int n_points,
693 const unsigned int base_entry = index * n_points;
695 static const std::array<unsigned int, 5> table = {{0, 1, 3, 2, 4}};
697 for (
unsigned int i = 0; i < n_points; ++i)
708 const unsigned int dimension,
709 const unsigned int set_num,
712 unsigned int new_dim;
731 for (
unsigned int d = 0; d < new_dim; ++d)
734 data_sets.back()[r * new_dim + d] = data_vectors(set_num + d, i);
750 const char *gmv_cell_type[4] = {
"",
"line 2",
"quad 4",
"hex 8"};
752 const char *ucd_cell_type[4] = {
"pt",
"line",
"quad",
"hex"};
754 const char *tecplot_cell_type[4] = {
"",
"lineseg",
"quadrilateral",
"brick"};
769 template <
int dim,
int spacedim>
770 std::array<unsigned int, 3>
772 const bool write_higher_order_cells)
774 std::array<unsigned int, 3> vtk_cell_id = {
780 if (write_higher_order_cells)
783 vtk_cell_id[2] = patch.
data.n_cols();
789 patch.
data.n_cols() == 6)
793 vtk_cell_id[2] = patch.
data.n_cols();
796 patch.
data.n_cols() == 10)
800 vtk_cell_id[2] = patch.
data.n_cols();
811 vtk_cell_id[2] = patch.
data.n_cols();
830 template <
int dim,
int spacedim>
832 get_equispaced_location(
834 const std::initializer_list<unsigned int> &lattice_location,
835 const unsigned int n_subdivisions)
842 const unsigned int xstep = (dim > 0 ? *(lattice_location.begin() + 0) : 0);
843 const unsigned int ystep = (dim > 1 ? *(lattice_location.begin() + 1) : 0);
844 const unsigned int zstep = (dim > 2 ? *(lattice_location.begin() + 2) : 0);
852 unsigned int point_no = 0;
857 point_no += (n_subdivisions + 1) * (n_subdivisions + 1) * zstep;
861 point_no += (n_subdivisions + 1) * ystep;
875 for (
unsigned int d = 0;
d < spacedim; ++
d)
876 node[d] = patch.
data(patch.
data.size(0) - spacedim + d, point_no);
883 if constexpr (dim == 0)
888 const double stepsize = 1. / n_subdivisions;
889 const double xfrac = xstep * stepsize;
895 const double yfrac = ystep * stepsize;
897 node += ((patch.
vertices[3] * xfrac) +
898 (patch.
vertices[2] * (1 - xfrac))) *
902 const double zfrac = zstep * stepsize;
904 node += (((patch.
vertices[5] * xfrac) +
905 (patch.
vertices[4] * (1 - xfrac))) *
908 (patch.
vertices[6] * (1 - xfrac))) *
921 template <
int dim,
int spacedim>
924 const unsigned int node_index)
929 unsigned int point_no_actual = node_index;
934 static const std::array<unsigned int, 5> table = {{0, 1, 3, 2, 4}};
935 point_no_actual = table[node_index];
943 for (
unsigned int d = 0;
d < spacedim; ++
d)
945 patch.
data(patch.
data.size(0) - spacedim + d, point_no_actual);
958 return patch.
vertices[point_no_actual];
969 template <
int dim,
int spacedim>
970 std::tuple<unsigned int, unsigned int>
971 count_nodes_and_cells(
974 unsigned int n_nodes = 0;
976 for (
const auto &patch : patches)
980 "The reference cell for this patch is set to 'Invalid', "
981 "but that is clearly not a valid choice. Did you forget "
982 "to set the reference cell for the patch?"));
997 return std::make_tuple(n_nodes, n_cells);
1007 template <
int dim,
int spacedim>
1008 std::tuple<unsigned int, unsigned int, unsigned int>
1009 count_nodes_and_cells_and_points(
1011 const bool write_higher_order_cells)
1013 unsigned int n_nodes = 0;
1015 unsigned int n_points_and_n_cells = 0;
1017 for (
const auto &patch : patches)
1021 n_nodes += Utilities::fixed_power<dim>(patch.
n_subdivisions + 1);
1023 if (write_higher_order_cells)
1029 n_points_and_n_cells +=
1038 const unsigned int n_subcells =
1041 n_points_and_n_cells +=
1047 n_nodes += patch.
data.n_cols();
1049 n_points_and_n_cells += patch.
data.n_cols() + 1;
1053 return std::make_tuple(n_nodes, n_cells, n_points_and_n_cells);
1061 template <
typename FlagsType>
1068 StreamBase(std::ostream &stream,
const FlagsType &flags)
1080 write_point(
const unsigned int,
const Point<dim> &)
1083 ExcMessage(
"The derived class you are using needs to "
1084 "reimplement this function if you want to call "
1104 write_cell(
const unsigned int ,
1105 const unsigned int ,
1106 std::array<unsigned int, dim> & )
1109 ExcMessage(
"The derived class you are using needs to "
1110 "reimplement this function if you want to call "
1122 write_cell_single(
const unsigned int index,
1123 const unsigned int start,
1124 const unsigned int n_points,
1130 (void)reference_cell;
1133 ExcMessage(
"The derived class you are using needs to "
1134 "reimplement this function if you want to call "
1152 template <
typename T>
1166 unsigned int selected_component;
1173 std::ostream &stream;
1178 const FlagsType flags;
1184 class DXStream :
public StreamBase<DataOutBase::DXFlags>
1191 write_point(
const unsigned int index,
const Point<dim> &);
1203 write_cell(
const unsigned int index,
1204 const unsigned int start,
1205 const std::array<unsigned int, dim> &offsets);
1213 template <
typename data>
1215 write_dataset(
const unsigned int index,
const std::vector<data> &values);
1221 class GmvStream :
public StreamBase<DataOutBase::GmvFlags>
1228 write_point(
const unsigned int index,
const Point<dim> &);
1240 write_cell(
const unsigned int index,
1241 const unsigned int start,
1242 const std::array<unsigned int, dim> &offsets);
1248 class TecplotStream :
public StreamBase<DataOutBase::TecplotFlags>
1255 write_point(
const unsigned int index,
const Point<dim> &);
1267 write_cell(
const unsigned int index,
1268 const unsigned int start,
1269 const std::array<unsigned int, dim> &offsets);
1275 class UcdStream :
public StreamBase<DataOutBase::UcdFlags>
1282 write_point(
const unsigned int index,
const Point<dim> &);
1296 write_cell(
const unsigned int index,
1297 const unsigned int start,
1298 const std::array<unsigned int, dim> &offsets);
1306 template <
typename data>
1308 write_dataset(
const unsigned int index,
const std::vector<data> &values);
1314 class VtkStream :
public StreamBase<DataOutBase::VtkFlags>
1321 write_point(
const unsigned int index,
const Point<dim> &);
1333 write_cell(
const unsigned int index,
1334 const unsigned int start,
1335 const std::array<unsigned int, dim> &offsets);
1342 write_cell_single(
const unsigned int index,
1343 const unsigned int start,
1344 const unsigned int n_points,
1356 write_high_order_cell(
const unsigned int start,
1357 const std::vector<unsigned> &connectivity);
1370 DXStream::write_point(
const unsigned int,
const Point<dim> &p)
1372 if (flags.coordinates_binary)
1375 for (
unsigned int d = 0;
d < dim; ++
d)
1377 stream.write(
reinterpret_cast<const char *
>(
data), dim *
sizeof(*
data));
1381 for (
unsigned int d = 0;
d < dim; ++
d)
1382 stream << p[d] <<
'\t';
1396 std::array<unsigned int, GeometryInfo<0>::vertices_per_cell>
1397 set_node_numbers(
const unsigned int ,
1398 const std::array<unsigned int, 0> & )
1406 std::array<unsigned int, GeometryInfo<1>::vertices_per_cell>
1407 set_node_numbers(
const unsigned int start,
1408 const std::array<unsigned int, 1> &offsets)
1410 std::array<unsigned int, GeometryInfo<1>::vertices_per_cell> nodes;
1412 nodes[1] = start + offsets[0];
1418 std::array<unsigned int, GeometryInfo<2>::vertices_per_cell>
1419 set_node_numbers(
const unsigned int start,
1420 const std::array<unsigned int, 2> &offsets)
1423 const unsigned int d1 = offsets[0];
1424 const unsigned int d2 = offsets[1];
1426 std::array<unsigned int, GeometryInfo<2>::vertices_per_cell> nodes;
1428 nodes[1] = start + d1;
1429 nodes[2] = start + d2;
1430 nodes[3] = start + d2 + d1;
1436 std::array<unsigned int, GeometryInfo<3>::vertices_per_cell>
1437 set_node_numbers(
const unsigned int start,
1438 const std::array<unsigned int, 3> &offsets)
1440 const unsigned int d1 = offsets[0];
1441 const unsigned int d2 = offsets[1];
1442 const unsigned int d3 = offsets[2];
1444 std::array<unsigned int, GeometryInfo<3>::vertices_per_cell> nodes;
1446 nodes[1] = start + d1;
1447 nodes[2] = start + d2;
1448 nodes[3] = start + d2 + d1;
1449 nodes[4] = start + d3;
1450 nodes[5] = start + d3 + d1;
1451 nodes[6] = start + d3 + d2;
1452 nodes[7] = start + d3 + d2 + d1;
1461 DXStream::write_cell(
const unsigned int,
1462 const unsigned int start,
1463 const std::array<unsigned int, dim> &offsets)
1466 DataOutBaseImplementation::set_node_numbers(start, offsets);
1468 if (flags.int_binary)
1470 std::array<unsigned int, GeometryInfo<dim>::vertices_per_cell> temp;
1471 for (
unsigned int i = 0; i < nodes.size(); ++i)
1473 stream.write(
reinterpret_cast<const char *
>(temp.data()),
1474 temp.size() *
sizeof(temp[0]));
1478 for (
unsigned int i = 0; i < nodes.size() - 1; ++i)
1480 stream << nodes[GeometryInfo<dim>::dx_to_deal[nodes.size() - 1]]
1487 template <
typename data>
1489 DXStream::write_dataset(
const unsigned int,
const std::vector<data> &values)
1491 if (flags.data_binary)
1493 stream.write(
reinterpret_cast<const char *
>(
values.data()),
1498 for (
unsigned int i = 0; i <
values.size(); ++i)
1499 stream <<
'\t' << values[i];
1515 GmvStream::write_point(
const unsigned int,
const Point<dim> &p)
1519 stream << p[selected_component] <<
' ';
1526 GmvStream::write_cell(
const unsigned int,
1527 const unsigned int s,
1528 const std::array<unsigned int, dim> &offsets)
1531 const unsigned int start = s + 1;
1532 stream << gmv_cell_type[dim] <<
'\n';
1544 const unsigned int d1 = offsets[0];
1546 stream <<
'\t' << start + d1;
1552 const unsigned int d1 = offsets[0];
1553 const unsigned int d2 = offsets[1];
1555 stream <<
'\t' << start + d1;
1556 stream <<
'\t' << start + d2 + d1 <<
'\t' << start + d2;
1562 const unsigned int d1 = offsets[0];
1563 const unsigned int d2 = offsets[1];
1564 const unsigned int d3 = offsets[2];
1566 stream <<
'\t' << start + d1;
1567 stream <<
'\t' << start + d2 + d1 <<
'\t' << start + d2;
1568 stream <<
'\t' << start + d3 <<
'\t' << start + d3 + d1 <<
'\t'
1569 << start + d3 + d2 + d1 <<
'\t' << start + d3 + d2;
1581 TecplotStream::TecplotStream(std::ostream &out,
1589 TecplotStream::write_point(
const unsigned int,
const Point<dim> &p)
1593 stream << p[selected_component] <<
'\n';
1600 TecplotStream::write_cell(
const unsigned int,
1601 const unsigned int s,
1602 const std::array<unsigned int, dim> &offsets)
1604 const unsigned int start = s + 1;
1616 const unsigned int d1 = offsets[0];
1618 stream <<
'\t' << start + d1;
1624 const unsigned int d1 = offsets[0];
1625 const unsigned int d2 = offsets[1];
1627 stream <<
'\t' << start + d1;
1628 stream <<
'\t' << start + d2 + d1 <<
'\t' << start + d2;
1634 const unsigned int d1 = offsets[0];
1635 const unsigned int d2 = offsets[1];
1636 const unsigned int d3 = offsets[2];
1638 stream <<
'\t' << start + d1;
1639 stream <<
'\t' << start + d2 + d1 <<
'\t' << start + d2;
1640 stream <<
'\t' << start + d3 <<
'\t' << start + d3 + d1 <<
'\t'
1641 << start + d3 + d2 + d1 <<
'\t' << start + d3 + d2;
1660 UcdStream::write_point(
const unsigned int index,
const Point<dim> &p)
1662 stream <<
index + 1 <<
" ";
1664 for (
unsigned int i = 0; i < dim; ++i)
1665 stream << p[i] <<
' ';
1667 for (
unsigned int i = dim; i < 3; ++i)
1676 UcdStream::write_cell(
const unsigned int index,
1677 const unsigned int start,
1678 const std::array<unsigned int, dim> &offsets)
1681 DataOutBaseImplementation::set_node_numbers(start, offsets);
1684 stream <<
index + 1 <<
"\t0 " << ucd_cell_type[dim];
1685 for (
unsigned int i = 0; i < nodes.size(); ++i)
1692 template <
typename data>
1694 UcdStream::write_dataset(
const unsigned int index,
1695 const std::vector<data> &values)
1697 stream <<
index + 1;
1698 for (
unsigned int i = 0; i <
values.size(); ++i)
1699 stream <<
'\t' << values[i];
1714 VtkStream::write_point(
const unsigned int,
const Point<dim> &p)
1719 for (
unsigned int i = dim; i < 3; ++i)
1728 VtkStream::write_cell(
const unsigned int,
1729 const unsigned int start,
1730 const std::array<unsigned int, dim> &offsets)
1732 stream << GeometryInfo<dim>::vertices_per_cell <<
'\t';
1744 const unsigned int d1 = offsets[0];
1746 stream <<
'\t' << start + d1;
1752 const unsigned int d1 = offsets[0];
1753 const unsigned int d2 = offsets[1];
1755 stream <<
'\t' << start + d1;
1756 stream <<
'\t' << start + d2 + d1 <<
'\t' << start + d2;
1762 const unsigned int d1 = offsets[0];
1763 const unsigned int d2 = offsets[1];
1764 const unsigned int d3 = offsets[2];
1766 stream <<
'\t' << start + d1;
1767 stream <<
'\t' << start + d2 + d1 <<
'\t' << start + d2;
1768 stream <<
'\t' << start + d3 <<
'\t' << start + d3 + d1 <<
'\t'
1769 << start + d3 + d2 + d1 <<
'\t' << start + d3 + d2;
1783 VtkStream::write_cell_single(
const unsigned int index,
1784 const unsigned int start,
1785 const unsigned int n_points,
1790 static const std::array<unsigned int, 5> table = {{0, 1, 3, 2, 4}};
1792 stream <<
'\t' << n_points;
1793 for (
unsigned int i = 0; i < n_points; ++i)
1804 VtkStream::write_high_order_cell(
const unsigned int start,
1805 const std::vector<unsigned> &connectivity)
1807 stream << connectivity.size();
1808 for (
const auto &c : connectivity)
1809 stream <<
'\t' << start + c;
1821 template <
int dim,
int spacedim>
1825 template <
int dim,
int spacedim>
1829 template <
int dim,
int spacedim>
1831 : patch_index(no_neighbor)
1833 , points_are_available(false)
1847 template <
int dim,
int spacedim>
1855 const double epsilon = 3e-16;
1857 if (vertices[i].distance(patch.
vertices[i]) > epsilon)
1873 if (
data.n_rows() != patch.
data.n_rows())
1876 if (
data.n_cols() != patch.
data.n_cols())
1879 for (
unsigned int i = 0; i <
data.n_rows(); ++i)
1880 for (
unsigned int j = 0; j <
data.n_cols(); ++j)
1881 if (
data[i][j] != patch.
data[i][j])
1889 template <
int dim,
int spacedim>
1893 return (
sizeof(vertices) /
sizeof(vertices[0]) *
1895 sizeof(neighbors) /
sizeof(neighbors[0]) *
1901 sizeof(reference_cell));
1906 template <
int dim,
int spacedim>
1910 std::swap(vertices, other_patch.vertices);
1911 std::swap(neighbors, other_patch.neighbors);
1912 std::swap(patch_index, other_patch.patch_index);
1913 std::swap(n_subdivisions, other_patch.n_subdivisions);
1914 data.swap(other_patch.data);
1915 std::swap(points_are_available, other_patch.points_are_available);
1916 std::swap(reference_cell, other_patch.reference_cell);
1921 template <
int spacedim>
1925 template <
int spacedim>
1929 template <
int spacedim>
1933 template <
int spacedim>
1936 template <
int spacedim>
1940 template <
int spacedim>
1942 : patch_index(no_neighbor)
1943 , points_are_available(false)
1950 template <
int spacedim>
1954 const unsigned int dim = 0;
1957 const double epsilon = 3e-16;
1959 if (vertices[i].distance(patch.
vertices[i]) > epsilon)
1968 if (
data.n_rows() != patch.
data.n_rows())
1971 if (
data.n_cols() != patch.
data.n_cols())
1974 for (
unsigned int i = 0; i <
data.n_rows(); ++i)
1975 for (
unsigned int j = 0; j <
data.n_cols(); ++j)
1976 if (
data[i][j] != patch.
data[i][j])
1984 template <
int spacedim>
1988 return (
sizeof(vertices) /
sizeof(vertices[0]) *
1996 template <
int spacedim>
2000 std::swap(vertices, other_patch.vertices);
2001 std::swap(patch_index, other_patch.patch_index);
2002 data.swap(other_patch.data);
2003 std::swap(points_are_available, other_patch.points_are_available);
2009 : write_preamble(write_preamble)
2024 : space_dimension_labels(labels)
2038 const bool bicubic_patch,
2039 const bool external_data)
2041 , bicubic_patch(bicubic_patch)
2042 , external_data(external_data)
2047 const bool xdmf_hdf5_output)
2048 : filter_duplicate_vertices(filter_duplicate_vertices)
2049 , xdmf_hdf5_output(xdmf_hdf5_output)
2057 "Filter duplicate vertices",
2060 "Whether to remove duplicate vertex values. deal.II duplicates "
2061 "vertices once for each adjacent cell so that it can output "
2062 "discontinuous quantities for which there may be more than one "
2063 "value for each vertex position. Setting this flag to "
2064 "'true' will merge all of these values by selecting a "
2065 "random one and outputting this as 'the' value for the vertex. "
2066 "As long as the data to be output corresponds to continuous "
2067 "fields, merging vertices has no effect. On the other hand, "
2068 "if the data to be output corresponds to discontinuous fields "
2069 "(either because you are using a discontinuous finite element, "
2070 "or because you are using a DataPostprocessor that yields "
2071 "discontinuous data, or because the data to be output has been "
2072 "produced by entirely different means), then the data in the "
2073 "output file no longer faithfully represents the underlying data "
2074 "because the discontinuous field has been replaced by a "
2075 "continuous one. Note also that the filtering can not occur "
2076 "on processor boundaries. Thus, a filtered discontinuous field "
2077 "looks like a continuous field inside of a subdomain, "
2078 "but like a discontinuous field at the subdomain boundary."
2080 "In any case, filtering results in drastically smaller output "
2081 "files (smaller by about a factor of 2^dim).");
2086 "Whether the data will be used in an XDMF/HDF5 combination.");
2101 const bool int_binary,
2102 const bool coordinates_binary,
2103 const bool data_binary)
2104 : write_neighbors(write_neighbors)
2105 , int_binary(int_binary)
2106 , coordinates_binary(coordinates_binary)
2107 , data_binary(data_binary)
2108 , data_double(false)
2118 "A boolean field indicating whether neighborship "
2119 "information between cells is to be written to the "
2120 "OpenDX output file");
2124 "Output format of integer numbers, which is "
2125 "either a text representation (ascii) or binary integer "
2126 "values of 32 or 64 bits length");
2130 "Output format of vertex coordinates, which is "
2131 "either a text representation (ascii) or binary "
2132 "floating point values of 32 or 64 bits length");
2136 "Output format of data values, which is "
2137 "either a text representation (ascii) or binary "
2138 "floating point values of 32 or 64 bits length");
2158 "A flag indicating whether a comment should be "
2159 "written to the beginning of the output file "
2160 "indicating date and time of creation as well "
2161 "as the creating program");
2175 const int azimuth_angle,
2176 const int polar_angle,
2177 const unsigned int line_thickness,
2179 const bool draw_colorbar)
2182 , height_vector(height_vector)
2183 , azimuth_angle(azimuth_angle)
2184 , polar_angle(polar_angle)
2185 , line_thickness(line_thickness)
2187 , draw_colorbar(draw_colorbar)
2198 "A flag indicating whether POVRAY should use smoothed "
2199 "triangles instead of the usual ones");
2203 "Whether POVRAY should use bicubic patches");
2207 "Whether camera and lighting information should "
2208 "be put into an external file \"data.inc\" or into "
2209 "the POVRAY input file");
2225 const unsigned int color_vector,
2227 const unsigned int size,
2228 const double line_width,
2229 const double azimut_angle,
2230 const double turn_angle,
2231 const double z_scaling,
2232 const bool draw_mesh,
2233 const bool draw_cells,
2234 const bool shade_cells,
2236 : height_vector(height_vector)
2237 , color_vector(color_vector)
2238 , size_type(size_type)
2240 , line_width(line_width)
2241 , azimut_angle(azimut_angle)
2242 , turn_angle(turn_angle)
2243 , z_scaling(z_scaling)
2244 , draw_mesh(draw_mesh)
2245 , draw_cells(draw_cells)
2246 , shade_cells(shade_cells)
2247 , color_function(color_function)
2286 double sum = xmax + xmin;
2287 double sum13 = xmin + 3 * xmax;
2288 double sum22 = 2 * xmin + 2 * xmax;
2289 double sum31 = 3 * xmin + xmax;
2290 double dif = xmax - xmin;
2291 double rezdif = 1.0 / dif;
2295 if (x < (sum31) / 4)
2297 else if (x < (sum22) / 4)
2299 else if (x < (sum13) / 4)
2310 rgb_values.
green = 0;
2311 rgb_values.
blue = (x - xmin) * 4. * rezdif;
2315 rgb_values.
green = (4 * x - 3 * xmin - xmax) * rezdif;
2316 rgb_values.
blue = (sum22 - 4. * x) * rezdif;
2319 rgb_values.
red = (4 * x - 2 * sum) * rezdif;
2320 rgb_values.
green = (xmin + 3 * xmax - 4 * x) * rezdif;
2321 rgb_values.
blue = 0;
2325 rgb_values.
green = (4 * x - xmin - 3 * xmax) * rezdif;
2326 rgb_values.
blue = (4. * x - sum13) * rezdif;
2333 rgb_values.
red = rgb_values.
green = rgb_values.
blue = 1;
2347 (x - xmin) / (xmax - xmin);
2360 1 - (x - xmin) / (xmax - xmin);
2372 "Number of the input vector that is to be used to "
2373 "generate height information");
2377 "Number of the input vector that is to be used to "
2378 "generate color information");
2382 "Whether width or height should be scaled to match "
2387 "The size (width or height) to which the eps output "
2388 "file is to be scaled");
2392 "The width in which the postscript renderer is to "
2397 "Angle of the viewing position against the vertical "
2402 "Angle of the viewing direction against the y-axis");
2406 "Scaling for the z-direction relative to the scaling "
2407 "used in x- and y-directions");
2411 "Whether the mesh lines, or only the surface should be "
2416 "Whether only the mesh lines, or also the interior of "
2417 "cells should be plotted. If this flag is false, then "
2418 "one can see through the mesh");
2422 "Whether the interior of cells shall be shaded");
2426 "default|grey scale|reverse grey scale"),
2427 "Name of a color function used to colorize mesh lines "
2428 "and/or cell interiors");
2438 if (prm.
get(
"Scale to width or height") ==
"width")
2450 if (prm.
get(
"Color function") ==
"default")
2452 else if (prm.
get(
"Color function") ==
"grey scale")
2454 else if (prm.
get(
"Color function") ==
"reverse grey scale")
2464 const bool keep_existing_file,
2465 const std::map<std::string,
2466 std::vector<std::pair<std::string, AttributeValue>>>
2469 , keep_existing_file(keep_existing_file)
2470 , attributes(attributes)
2474 : compression_level(compression_level)
2479 : zone_name(zone_name)
2480 , solution_time(solution_time)
2494 const unsigned int cycle,
2495 const bool print_date_and_time,
2497 const bool write_higher_order_cells,
2498 const std::map<std::string, std::string> &physical_units)
2501 , print_date_and_time(print_date_and_time)
2502 , compression_level(compression_level)
2503 , write_higher_order_cells(write_higher_order_cells)
2504 , physical_units(physical_units)
2512 if (format_name ==
"none")
2515 if (format_name ==
"dx")
2518 if (format_name ==
"ucd")
2521 if (format_name ==
"gnuplot")
2524 if (format_name ==
"povray")
2527 if (format_name ==
"eps")
2530 if (format_name ==
"gmv")
2533 if (format_name ==
"tecplot")
2536 if (format_name ==
"vtk")
2539 if (format_name ==
"vtu")
2542 if (format_name ==
"deal.II intermediate")
2545 if (format_name ==
"hdf5")
2549 ExcMessage(
"The given file format name is not recognized: <" +
2550 format_name +
">"));
2561 return "none|dx|ucd|gnuplot|povray|eps|gmv|tecplot|vtk|vtu|hdf5|svg|deal.II intermediate";
2569 switch (output_format)
2611 template <
int dim,
int spacedim>
2612 std::vector<Point<spacedim>>
2616 static const std::array<unsigned int, 5> table = {{0, 1, 3, 2, 4}};
2618 std::vector<Point<spacedim>> node_positions;
2619 std::size_t n_nodes = 0;
2620 for (
const auto &patch : patches)
2622 if (patch.
reference_cell != ReferenceCells::get_hypercube<dim>())
2623 n_nodes += patch.
data.n_cols();
2625 n_nodes += Utilities::fixed_power<dim>(patch.
n_subdivisions + 1);
2627 node_positions.reserve(n_nodes);
2629 for (
const auto &patch : patches)
2632 if (patch.
reference_cell != ReferenceCells::get_hypercube<dim>())
2634 for (
unsigned int point_no = 0; point_no < patch.
data.n_cols();
2636 node_positions.emplace_back(get_node_location(
2645 const unsigned int n = n_subdivisions + 1;
2650 node_positions.emplace_back(
2651 get_equispaced_location(patch, {}, n_subdivisions));
2654 for (
unsigned int i1 = 0; i1 < n; ++i1)
2655 node_positions.emplace_back(
2656 get_equispaced_location(patch, {i1}, n_subdivisions));
2659 for (
unsigned int i2 = 0; i2 < n; ++i2)
2660 for (
unsigned int i1 = 0; i1 < n; ++i1)
2661 node_positions.emplace_back(get_equispaced_location(
2662 patch, {i1, i2}, n_subdivisions));
2665 for (
unsigned int i3 = 0; i3 < n; ++i3)
2666 for (
unsigned int i2 = 0; i2 < n; ++i2)
2667 for (
unsigned int i1 = 0; i1 < n; ++i1)
2668 node_positions.emplace_back(get_equispaced_location(
2669 patch, {i1, i2, i3}, n_subdivisions));
2678 return node_positions;
2682 template <
int dim,
int spacedim,
typename StreamType>
2688 const std::vector<Point<spacedim>> node_positions =
2692 for (
const auto &node : node_positions)
2693 out.write_point(count++, node);
2699 template <
int dim,
int spacedim,
typename StreamType>
2704 unsigned int count = 0;
2705 unsigned int first_vertex_of_patch = 0;
2706 for (
const auto &patch : patches)
2709 if (patch.
reference_cell != ReferenceCells::get_hypercube<dim>())
2711 out.write_cell_single(count++,
2712 first_vertex_of_patch,
2713 patch.
data.n_cols(),
2715 first_vertex_of_patch += patch.
data.n_cols();
2720 const unsigned int n = n_subdivisions + 1;
2726 const unsigned int offset = first_vertex_of_patch;
2727 out.template write_cell<0>(count++, offset, {});
2733 constexpr unsigned int d1 = 1;
2735 for (
unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
2737 const unsigned int offset =
2738 first_vertex_of_patch + i1 * d1;
2739 out.template write_cell<1>(count++, offset, {{d1}});
2747 constexpr unsigned int d1 = 1;
2748 const unsigned int d2 = n;
2750 for (
unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
2751 for (
unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
2753 const unsigned int offset =
2754 first_vertex_of_patch + i2 * d2 + i1 * d1;
2755 out.template write_cell<2>(count++,
2765 constexpr unsigned int d1 = 1;
2766 const unsigned int d2 = n;
2767 const unsigned int d3 = n * n;
2769 for (
unsigned int i3 = 0; i3 < n_subdivisions; ++i3)
2770 for (
unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
2771 for (
unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
2773 const unsigned int offset = first_vertex_of_patch +
2776 out.template write_cell<3>(count++,
2788 first_vertex_of_patch +=
2789 Utilities::fixed_power<dim>(n_subdivisions + 1);
2798 template <
int dim,
int spacedim,
typename StreamType>
2802 const bool legacy_format)
2806 (void)legacy_format;
2808 unsigned int first_vertex_of_patch = 0;
2810 std::vector<unsigned> connectivity;
2812 for (
const auto &patch : patches)
2814 if (patch.
reference_cell != ReferenceCells::get_hypercube<dim>())
2816 connectivity.resize(patch.
data.n_cols());
2818 for (
unsigned int i = 0; i < patch.
data.n_cols(); ++i)
2819 connectivity[i] = i;
2821 out.template write_high_order_cell<dim>(first_vertex_of_patch,
2824 first_vertex_of_patch += patch.
data.n_cols();
2829 const unsigned int n = n_subdivisions + 1;
2831 connectivity.resize(Utilities::fixed_power<dim>(n));
2833 if constexpr (dim == 0)
2836 ExcMessage(
"Point-like cells should not be possible "
2837 "when writing higher-order cells."));
2839 else if constexpr (dim == 1)
2841 for (
unsigned int i1 = 0; i1 < n_subdivisions + 1; ++i1)
2843 const unsigned int local_index = i1;
2844 const unsigned int connectivity_index =
2846 {{i1}}, {{n_subdivisions}}, legacy_format);
2847 connectivity[connectivity_index] = local_index;
2850 else if constexpr (dim == 2)
2852 for (
unsigned int i2 = 0; i2 < n_subdivisions + 1; ++i2)
2853 for (
unsigned int i1 = 0; i1 < n_subdivisions + 1; ++i1)
2855 const unsigned int local_index = i2 * n + i1;
2856 const unsigned int connectivity_index =
2859 {{n_subdivisions, n_subdivisions}},
2861 connectivity[connectivity_index] = local_index;
2864 else if constexpr (dim == 3)
2866 for (
unsigned int i3 = 0; i3 < n_subdivisions + 1; ++i3)
2867 for (
unsigned int i2 = 0; i2 < n_subdivisions + 1; ++i2)
2868 for (
unsigned int i1 = 0; i1 < n_subdivisions + 1; ++i1)
2870 const unsigned int local_index =
2871 i3 * n * n + i2 * n + i1;
2872 const unsigned int connectivity_index =
2875 {{n_subdivisions, n_subdivisions, n_subdivisions}},
2877 connectivity[connectivity_index] = local_index;
2885 out.template write_high_order_cell<dim>(first_vertex_of_patch,
2889 first_vertex_of_patch += Utilities::fixed_power<dim>(n);
2897 template <
int dim,
int spacedim,
typename StreamType>
2900 unsigned int n_data_sets,
2901 const bool double_precision,
2905 unsigned int count = 0;
2907 for (
const auto &patch : patches)
2910 const unsigned int n = n_subdivisions + 1;
2912 Assert((patch.
data.n_rows() == n_data_sets &&
2914 (patch.
data.n_rows() == n_data_sets + spacedim &&
2917 (n_data_sets + spacedim) :
2919 patch.
data.n_rows()));
2920 Assert(patch.
data.n_cols() == Utilities::fixed_power<dim>(n),
2923 std::vector<float> floats(n_data_sets);
2924 std::vector<double> doubles(n_data_sets);
2927 for (
unsigned int i = 0; i < Utilities::fixed_power<dim>(n);
2929 if (double_precision)
2931 for (
unsigned int data_set = 0; data_set < n_data_sets;
2933 doubles[data_set] = patch.
data(data_set, i);
2934 out.write_dataset(count, doubles);
2938 for (
unsigned int data_set = 0; data_set < n_data_sets;
2940 floats[data_set] = patch.
data(data_set, i);
2941 out.write_dataset(count, floats);
2966 camera_vertical[0] = camera_horizontal[1] * camera_direction[2] -
2967 camera_horizontal[2] * camera_direction[1];
2968 camera_vertical[1] = camera_horizontal[2] * camera_direction[0] -
2969 camera_horizontal[0] * camera_direction[2];
2970 camera_vertical[2] = camera_horizontal[0] * camera_direction[1] -
2971 camera_horizontal[1] * camera_direction[0];
2975 phi /= (point[0] - camera_position[0]) * camera_direction[0] +
2976 (point[1] - camera_position[1]) * camera_direction[1] +
2977 (point[2] - camera_position[2]) * camera_direction[2];
2981 camera_position[0] + phi * (point[0] - camera_position[0]);
2983 camera_position[1] + phi * (point[1] - camera_position[1]);
2985 camera_position[2] + phi * (point[2] - camera_position[2]);
2988 projection_decomposition[0] = (projection[0] - camera_position[0] -
2989 camera_focus * camera_direction[0]) *
2990 camera_horizontal[0];
2991 projection_decomposition[0] += (projection[1] - camera_position[1] -
2992 camera_focus * camera_direction[1]) *
2993 camera_horizontal[1];
2994 projection_decomposition[0] += (projection[2] - camera_position[2] -
2995 camera_focus * camera_direction[2]) *
2996 camera_horizontal[2];
2998 projection_decomposition[1] = (projection[0] - camera_position[0] -
2999 camera_focus * camera_direction[0]) *
3001 projection_decomposition[1] += (projection[1] - camera_position[1] -
3002 camera_focus * camera_direction[1]) *
3004 projection_decomposition[1] += (projection[2] - camera_position[2] -
3005 camera_focus * camera_direction[2]) *
3008 return projection_decomposition;
3017 svg_get_gradient_parameters(
Point<3> points[])
3023 for (
int i = 0; i < 2; ++i)
3025 for (
int j = 0; j < 2 - i; ++j)
3027 if (points[j][2] > points[j + 1][2])
3030 points[j] = points[j + 1];
3031 points[j + 1] = temp;
3038 v_inter = points[1];
3045 A[0][0] = v_max[0] - v_min[0];
3046 A[0][1] = v_inter[0] - v_min[0];
3047 A[1][0] = v_max[1] - v_min[1];
3048 A[1][1] = v_inter[1] - v_min[1];
3054 bool col_change =
false;
3063 double temp = A[1][0];
3068 for (
unsigned int k = 0; k < 1; ++k)
3070 for (
unsigned int i = k + 1; i < 2; ++i)
3072 x = A[i][k] / A[k][k];
3074 for (
unsigned int j = k + 1; j < 2; ++j)
3075 A[i][j] = A[i][j] - A[k][j] * x;
3077 b[i] = b[i] - b[k] * x;
3081 b[1] =
b[1] /
A[1][1];
3083 for (
int i = 0; i >= 0; i--)
3087 for (
unsigned int j = i + 1; j < 2; ++j)
3090 b[i] =
sum /
A[i][i];
3100 double c =
b[0] * (v_max[2] - v_min[2]) + b[1] * (v_inter[2] - v_min[2]) +
3104 A[0][0] = v_max[0] - v_min[0];
3105 A[0][1] = v_inter[0] - v_min[0];
3106 A[1][0] = v_max[1] - v_min[1];
3107 A[1][1] = v_inter[1] - v_min[1];
3109 b[0] = 1.0 - v_min[0];
3121 double temp = A[1][0];
3126 for (
unsigned int k = 0; k < 1; ++k)
3128 for (
unsigned int i = k + 1; i < 2; ++i)
3130 x = A[i][k] / A[k][k];
3132 for (
unsigned int j = k + 1; j < 2; ++j)
3133 A[i][j] = A[i][j] - A[k][j] * x;
3135 b[i] = b[i] - b[k] * x;
3139 b[1] =
b[1] /
A[1][1];
3141 for (
int i = 0; i >= 0; i--)
3145 for (
unsigned int j = i + 1; j < 2; ++j)
3146 sum = sum - A[i][j] * b[j];
3148 b[i] = sum / A[i][i];
3158 gradient[0] = b[0] * (v_max[2] - v_min[2]) +
3159 b[1] * (v_inter[2] - v_min[2]) - c + v_min[2];
3162 A[0][0] = v_max[0] - v_min[0];
3163 A[0][1] = v_inter[0] - v_min[0];
3164 A[1][0] = v_max[1] - v_min[1];
3165 A[1][1] = v_inter[1] - v_min[1];
3168 b[1] = 1.0 - v_min[1];
3179 double temp = A[1][0];
3184 for (
unsigned int k = 0; k < 1; ++k)
3186 for (
unsigned int i = k + 1; i < 2; ++i)
3188 x = A[i][k] / A[k][k];
3190 for (
unsigned int j = k + 1; j < 2; ++j)
3191 A[i][j] = A[i][j] - A[k][j] * x;
3193 b[i] = b[i] - b[k] * x;
3197 b[1] =
b[1] /
A[1][1];
3199 for (
int i = 0; i >= 0; i--)
3203 for (
unsigned int j = i + 1; j < 2; ++j)
3204 sum = sum - A[i][j] * b[j];
3206 b[i] = sum / A[i][i];
3216 gradient[1] = b[0] * (v_max[2] - v_min[2]) +
3217 b[1] * (v_inter[2] - v_min[2]) - c + v_min[2];
3220 gradient /= gradient.norm();
3222 const double lambda = -gradient[0] * (v_min[0] - v_max[0]) -
3223 gradient[1] * (v_min[1] - v_max[1]);
3227 gradient_parameters[0] = v_min[0];
3228 gradient_parameters[1] = v_min[1];
3230 gradient_parameters[2] = v_min[0] + lambda * gradient[0];
3231 gradient_parameters[3] = v_min[1] + lambda * gradient[1];
3233 gradient_parameters[4] = v_min[2];
3234 gradient_parameters[5] = v_max[2];
3236 return gradient_parameters;
3242 template <
int dim,
int spacedim>
3246 const std::vector<std::string> &data_names,
3248 std::tuple<
unsigned int,
3252 const UcdFlags &flags,
3261#ifndef DEAL_II_WITH_MPI
3270 if (patches.empty())
3274 const unsigned int n_data_sets = data_names.size();
3276 UcdStream ucd_out(out, flags);
3280 auto [n_nodes, n_cells] = count_nodes_and_cells(patches);
3283 if (flags.write_preamble)
3286 <<
"# This file was generated by the deal.II library." <<
'\n'
3290 <<
"# For a description of the UCD format see the AVS Developer's guide."
3296 out << n_nodes <<
' ' << n_cells <<
' ' << n_data_sets <<
' ' << 0
3309 if (n_data_sets != 0)
3311 out << n_data_sets <<
" ";
3312 for (
unsigned int i = 0; i < n_data_sets; ++i)
3317 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
3318 out << data_names[data_set]
3322 write_data(patches, n_data_sets,
true, ucd_out);
3332 template <
int dim,
int spacedim>
3336 const std::vector<std::string> &data_names,
3338 std::tuple<
unsigned int,
3350#ifndef DEAL_II_WITH_MPI
3359 if (patches.empty())
3363 DXStream dx_out(out, flags);
3366 unsigned int offset = 0;
3368 const unsigned int n_data_sets = data_names.size();
3372 auto [n_nodes, n_cells] = count_nodes_and_cells(patches);
3375 out <<
"object \"vertices\" class array type float rank 1 shape "
3376 << spacedim <<
" items " << n_nodes;
3380 out <<
" lsb ieee data 0" <<
'\n';
3381 offset += n_nodes * spacedim *
sizeof(float);
3385 out <<
" data follows" <<
'\n';
3394 out <<
"object \"cells\" class array type int rank 1 shape "
3399 out <<
" lsb binary data " << offset <<
'\n';
3400 offset += n_cells *
sizeof(
int);
3404 out <<
" data follows" <<
'\n';
3410 out <<
"attribute \"element type\" string \"";
3411 if constexpr (dim == 1)
3413 else if constexpr (dim == 2)
3415 else if constexpr (dim == 3)
3417 out <<
"\"" <<
'\n' <<
"attribute \"ref\" string \"positions\"" <<
'\n';
3424 out <<
"object \"neighbors\" class array type int rank 1 shape "
3428 for (
const auto &patch : patches)
3431 const unsigned int n1 = (dim > 0) ? n : 1;
3432 const unsigned int n2 = (dim > 1) ? n : 1;
3433 const unsigned int n3 = (dim > 2) ? n : 1;
3434 const unsigned int x_minus = (dim > 0) ? 0 : 0;
3435 const unsigned int x_plus = (dim > 0) ? 1 : 0;
3436 const unsigned int y_minus = (dim > 1) ? 2 : 0;
3437 const unsigned int y_plus = (dim > 1) ? 3 : 0;
3438 const unsigned int z_minus = (dim > 2) ? 4 : 0;
3439 const unsigned int z_plus = (dim > 2) ? 5 : 0;
3440 unsigned int cells_per_patch = Utilities::fixed_power<dim>(n);
3441 unsigned int dx = 1;
3442 unsigned int dy = n;
3443 unsigned int dz = n * n;
3445 const unsigned int patch_start =
3448 for (
unsigned int i3 = 0; i3 < n3; ++i3)
3449 for (
unsigned int i2 = 0; i2 < n2; ++i2)
3450 for (
unsigned int i1 = 0; i1 < n1; ++i1)
3452 const unsigned int nx = i1 *
dx;
3453 const unsigned int ny = i2 * dy;
3454 const unsigned int nz = i3 * dz;
3466 const unsigned int nn = patch.
neighbors[x_minus];
3470 << (nn * cells_per_patch + ny + nz +
dx * (n - 1));
3476 out <<
'\t' << patch_start + nx -
dx + ny + nz;
3481 const unsigned int nn = patch.
neighbors[x_plus];
3484 out << (nn * cells_per_patch + ny + nz);
3490 out <<
'\t' << patch_start + nx +
dx + ny + nz;
3497 const unsigned int nn = patch.
neighbors[y_minus];
3501 << (nn * cells_per_patch + nx + nz + dy * (n - 1));
3507 out <<
'\t' << patch_start + nx + ny - dy + nz;
3512 const unsigned int nn = patch.
neighbors[y_plus];
3515 out << (nn * cells_per_patch + nx + nz);
3521 out <<
'\t' << patch_start + nx + ny + dy + nz;
3529 const unsigned int nn = patch.
neighbors[z_minus];
3533 << (nn * cells_per_patch + nx + ny + dz * (n - 1));
3539 out <<
'\t' << patch_start + nx + ny + nz - dz;
3544 const unsigned int nn = patch.
neighbors[z_plus];
3547 out << (nn * cells_per_patch + nx + ny);
3553 out <<
'\t' << patch_start + nx + ny + nz + dz;
3561 if (n_data_sets != 0)
3563 out <<
"object \"data\" class array type float rank 1 shape "
3564 << n_data_sets <<
" items " << n_nodes;
3568 out <<
" lsb ieee data " << offset <<
'\n';
3569 offset += n_data_sets * n_nodes *
3570 ((flags.
data_double) ?
sizeof(
double) :
sizeof(float));
3574 out <<
" data follows" <<
'\n';
3579 out <<
"attribute \"dep\" string \"positions\"" <<
'\n';
3583 out <<
"object \"data\" class constantarray type float rank 0 items "
3584 << n_nodes <<
" data follows" <<
'\n'
3590 out <<
"object \"deal data\" class field" <<
'\n'
3591 <<
"component \"positions\" value \"vertices\"" <<
'\n'
3592 <<
"component \"connections\" value \"cells\"" <<
'\n'
3593 <<
"component \"data\" value \"data\"" <<
'\n';
3596 out <<
"component \"neighbors\" value \"neighbors\"" <<
'\n';
3603 out <<
"end" <<
'\n';
3621 template <
int dim,
int spacedim>
3625 const std::vector<std::string> &data_names,
3627 std::tuple<
unsigned int,
3636#ifndef DEAL_II_WITH_MPI
3646 if (patches.empty())
3650 const unsigned int n_data_sets = data_names.size();
3654 out <<
"# This file was generated by the deal.II library." <<
'\n'
3658 <<
"# For a description of the GNUPLOT format see the GNUPLOT manual."
3665 for (
unsigned int spacedim_n = 0; spacedim_n < spacedim; ++spacedim_n)
3670 for (
const auto &data_name : data_names)
3671 out <<
'<' << data_name <<
"> ";
3677 for (
const auto &patch : patches)
3680 const unsigned int n_points_per_direction = n_subdivisions + 1;
3682 Assert((patch.
data.n_rows() == n_data_sets &&
3684 (patch.
data.n_rows() == n_data_sets + spacedim &&
3687 (n_data_sets + spacedim) :
3689 patch.
data.n_rows()));
3691 auto output_point_data =
3692 [&out, &patch, n_data_sets](
const unsigned int point_index)
mutable {
3693 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
3694 out << patch.data(data_set, point_index) <<
' ';
3703 Assert(patch.data.n_cols() == 1,
3705 n_subdivisions + 1));
3709 out << get_equispaced_location(patch, {}, n_subdivisions)
3711 output_point_data(0);
3721 Assert(patch.data.n_cols() ==
3722 Utilities::fixed_power<dim>(n_points_per_direction),
3724 n_subdivisions + 1));
3726 for (
unsigned int i1 = 0; i1 < n_points_per_direction; ++i1)
3729 out << get_equispaced_location(patch, {i1}, n_subdivisions)
3732 output_point_data(i1);
3745 Assert(patch.data.n_cols() == Utilities::fixed_power<dim>(
3746 n_points_per_direction),
3748 n_subdivisions + 1));
3750 for (
unsigned int i2 = 0; i2 < n_points_per_direction; ++i2)
3752 for (
unsigned int i1 = 0; i1 < n_points_per_direction;
3756 out << get_equispaced_location(patch,
3761 output_point_data(i1 + i2 * n_points_per_direction);
3786 out << get_node_location(patch, 0) <<
' ';
3787 output_point_data(0);
3790 out << get_node_location(patch, 1) <<
' ';
3791 output_point_data(1);
3795 out << get_node_location(patch, 2) <<
' ';
3796 output_point_data(2);
3799 out << get_node_location(patch, 2) <<
' ';
3800 output_point_data(2);
3820 Assert(patch.data.n_cols() == Utilities::fixed_power<dim>(
3821 n_points_per_direction),
3823 n_subdivisions + 1));
3828 for (
unsigned int i3 = 0; i3 < n_points_per_direction; ++i3)
3829 for (
unsigned int i2 = 0; i2 < n_points_per_direction;
3831 for (
unsigned int i1 = 0; i1 < n_points_per_direction;
3836 get_equispaced_location(patch,
3840 if (i1 < n_subdivisions)
3843 out << this_point <<
' ';
3844 output_point_data(i1 +
3845 i2 * n_points_per_direction +
3846 i3 * n_points_per_direction *
3847 n_points_per_direction);
3851 out << get_equispaced_location(patch,
3856 output_point_data((i1 + 1) +
3857 i2 * n_points_per_direction +
3858 i3 * n_points_per_direction *
3859 n_points_per_direction);
3863 out <<
'\n' <<
'\n';
3867 if (i2 < n_subdivisions)
3870 out << this_point <<
' ';
3871 output_point_data(i1 +
3872 i2 * n_points_per_direction +
3873 i3 * n_points_per_direction *
3874 n_points_per_direction);
3878 out << get_equispaced_location(patch,
3884 i1 + (i2 + 1) * n_points_per_direction +
3885 i3 * n_points_per_direction *
3886 n_points_per_direction);
3890 out <<
'\n' <<
'\n';
3894 if (i3 < n_subdivisions)
3897 out << this_point <<
' ';
3898 output_point_data(i1 +
3899 i2 * n_points_per_direction +
3900 i3 * n_points_per_direction *
3901 n_points_per_direction);
3905 out << get_equispaced_location(patch,
3911 i1 + i2 * n_points_per_direction +
3912 (i3 + 1) * n_points_per_direction *
3913 n_points_per_direction);
3916 out <<
'\n' <<
'\n';
3925 for (
const unsigned int v : {0, 1, 2, 0, 3, 2})
3927 out << get_node_location(patch, v) <<
' ';
3928 output_point_data(v);
3933 for (
const unsigned int v : {3, 1})
3935 out << get_node_location(patch, v) <<
' ';
3936 output_point_data(v);
3946 for (
const unsigned int v : {0, 1, 3, 2, 0, 4, 1})
3948 out << get_node_location(patch, v) <<
' ';
3949 output_point_data(v);
3954 for (
const unsigned int v : {2, 4, 3})
3956 out << get_node_location(patch, v) <<
' ';
3957 output_point_data(v);
3971 for (
const unsigned int v : {0, 1, 2, 0, 3, 4, 5, 3})
3973 out << get_node_location(patch, v) <<
' ';
3974 output_point_data(v);
3979 for (
const unsigned int v : {1, 4})
3981 out << get_node_location(patch, v) <<
' ';
3982 output_point_data(v);
3987 for (
const unsigned int v : {2, 5})
3989 out << get_node_location(patch, v) <<
' ';
3990 output_point_data(v);
4015 template <
int dim,
int spacedim>
4017 do_write_povray(
const std::vector<Patch<dim, spacedim>> &,
4018 const std::vector<std::string> &,
4019 const PovrayFlags &,
4023 ExcMessage(
"Writing files in POVRAY format is only supported "
4024 "for two-dimensional meshes."));
4030 do_write_povray(
const std::vector<Patch<2, 2>> &patches,
4031 const std::vector<std::string> &data_names,
4032 const PovrayFlags &flags,
4037#ifndef DEAL_II_WITH_MPI
4046 if (patches.empty())
4049 constexpr int dim = 2;
4051 constexpr int spacedim = 2;
4053 const unsigned int n_data_sets = data_names.size();
4059 <<
"/* This file was generated by the deal.II library." <<
'\n'
4063 <<
" For a description of the POVRAY format see the POVRAY manual."
4068 out <<
"#include \"colors.inc\" " <<
'\n'
4069 <<
"#include \"textures.inc\" " <<
'\n';
4073 if (flags.external_data)
4074 out <<
"#include \"data.inc\" " <<
'\n';
4080 <<
"camera {" <<
'\n'
4081 <<
" location <1,4,-7>" <<
'\n'
4082 <<
" look_at <0,0,0>" <<
'\n'
4083 <<
" angle 30" <<
'\n'
4088 <<
"light_source {" <<
'\n'
4089 <<
" <1,4,-7>" <<
'\n'
4090 <<
" color Grey" <<
'\n'
4093 <<
"light_source {" <<
'\n'
4094 <<
" <0,20,0>" <<
'\n'
4095 <<
" color White" <<
'\n'
4102 double hmin = patches[0].data(0, 0);
4103 double hmax = patches[0].data(0, 0);
4105 for (
const auto &patch : patches)
4109 Assert((patch.
data.n_rows() == n_data_sets &&
4111 (patch.
data.n_rows() == n_data_sets + spacedim &&
4114 (n_data_sets + spacedim) :
4116 patch.
data.n_rows()));
4118 Utilities::fixed_power<dim>(n_subdivisions + 1),
4120 n_subdivisions + 1));
4122 for (
unsigned int i = 0; i < n_subdivisions + 1; ++i)
4123 for (
unsigned int j = 0; j < n_subdivisions + 1; ++j)
4125 const int dl = i * (n_subdivisions + 1) + j;
4126 if (patch.
data(0, dl) < hmin)
4127 hmin = patch.
data(0, dl);
4128 if (patch.
data(0, dl) > hmax)
4129 hmax = patch.
data(0, dl);
4133 out <<
"#declare HMIN=" << hmin <<
";" <<
'\n'
4134 <<
"#declare HMAX=" << hmax <<
";" <<
'\n'
4137 if (!flags.external_data)
4140 out <<
"#declare Tex=texture{" <<
'\n'
4141 <<
" pigment {" <<
'\n'
4142 <<
" gradient y" <<
'\n'
4143 <<
" scale y*(HMAX-HMIN)*" << 0.1 <<
'\n'
4144 <<
" color_map {" <<
'\n'
4145 <<
" [0.00 color Light_Purple] " <<
'\n'
4146 <<
" [0.95 color Light_Purple] " <<
'\n'
4147 <<
" [1.00 color White] " <<
'\n'
4152 if (!flags.bicubic_patch)
4155 out <<
'\n' <<
"mesh {" <<
'\n';
4159 for (
const auto &patch : patches)
4162 const unsigned int n = n_subdivisions + 1;
4163 const unsigned int d1 = 1;
4164 const unsigned int d2 = n;
4166 Assert((patch.
data.n_rows() == n_data_sets &&
4168 (patch.
data.n_rows() == n_data_sets + spacedim &&
4171 (n_data_sets + spacedim) :
4173 patch.
data.n_rows()));
4174 Assert(patch.
data.n_cols() == Utilities::fixed_power<dim>(n),
4176 n_subdivisions + 1));
4179 std::vector<Point<spacedim>> ver(n * n);
4181 for (
unsigned int i2 = 0; i2 < n; ++i2)
4182 for (
unsigned int i1 = 0; i1 < n; ++i1)
4185 ver[i1 * d1 + i2 * d2] =
4186 get_equispaced_location(patch, {i1, i2}, n_subdivisions);
4190 if (!flags.bicubic_patch)
4193 std::vector<Point<3>> nrml;
4203 for (
unsigned int i = 0; i < n; ++i)
4204 for (
unsigned int j = 0; j < n; ++j)
4206 const unsigned int il = (i == 0) ? i : (i - 1);
4207 const unsigned int ir =
4208 (i == n_subdivisions) ? i : (i + 1);
4209 const unsigned int jl = (j == 0) ? j : (j - 1);
4210 const unsigned int jr =
4211 (j == n_subdivisions) ? j : (j + 1);
4214 ver[ir * d1 + j * d2][0] - ver[il * d1 + j * d2][0];
4215 h1[1] = patch.
data(0, ir * d1 + j * d2) -
4216 patch.
data(0, il * d1 + j * d2);
4218 ver[ir * d1 + j * d2][1] - ver[il * d1 + j * d2][1];
4221 ver[i * d1 + jr * d2][0] - ver[i * d1 + jl * d2][0];
4222 h2[1] = patch.
data(0, i * d1 + jr * d2) -
4223 patch.
data(0, i * d1 + jl * d2);
4225 ver[i * d1 + jr * d2][1] - ver[i * d1 + jl * d2][1];
4227 nrml[i * d1 + j * d2][0] =
4228 h1[1] * h2[2] - h1[2] * h2[1];
4229 nrml[i * d1 + j * d2][1] =
4230 h1[2] * h2[0] - h1[0] * h2[2];
4231 nrml[i * d1 + j * d2][2] =
4232 h1[0] * h2[1] - h1[1] * h2[0];
4235 double norm = std::hypot(nrml[i * d1 + j * d2][0],
4236 nrml[i * d1 + j * d2][1],
4237 nrml[i * d1 + j * d2][2]);
4239 if (nrml[i * d1 + j * d2][1] < 0)
4242 for (
unsigned int k = 0; k < 3; ++k)
4243 nrml[i * d1 + j * d2][k] /= norm;
4248 for (
unsigned int i = 0; i < n_subdivisions; ++i)
4249 for (
unsigned int j = 0; j < n_subdivisions; ++j)
4252 const int dl = i * d1 + j * d2;
4258 out <<
"smooth_triangle {" <<
'\n'
4259 <<
"\t<" << ver[dl][0] <<
"," << patch.
data(0, dl)
4260 <<
"," << ver[dl][1] <<
">, <" << nrml[dl][0]
4261 <<
", " << nrml[dl][1] <<
", " << nrml[dl][2]
4263 out <<
" \t<" << ver[dl + d1][0] <<
","
4264 << patch.
data(0, dl + d1) <<
"," << ver[dl + d1][1]
4265 <<
">, <" << nrml[dl + d1][0] <<
", "
4266 << nrml[dl + d1][1] <<
", " << nrml[dl + d1][2]
4268 out <<
"\t<" << ver[dl + d1 + d2][0] <<
","
4269 << patch.
data(0, dl + d1 + d2) <<
","
4270 << ver[dl + d1 + d2][1] <<
">, <"
4271 << nrml[dl + d1 + d2][0] <<
", "
4272 << nrml[dl + d1 + d2][1] <<
", "
4273 << nrml[dl + d1 + d2][2] <<
">}" <<
'\n';
4276 out <<
"smooth_triangle {" <<
'\n'
4277 <<
"\t<" << ver[dl][0] <<
"," << patch.
data(0, dl)
4278 <<
"," << ver[dl][1] <<
">, <" << nrml[dl][0]
4279 <<
", " << nrml[dl][1] <<
", " << nrml[dl][2]
4281 out <<
"\t<" << ver[dl + d1 + d2][0] <<
","
4282 << patch.
data(0, dl + d1 + d2) <<
","
4283 << ver[dl + d1 + d2][1] <<
">, <"
4284 << nrml[dl + d1 + d2][0] <<
", "
4285 << nrml[dl + d1 + d2][1] <<
", "
4286 << nrml[dl + d1 + d2][2] <<
">," <<
'\n';
4287 out <<
"\t<" << ver[dl + d2][0] <<
","
4288 << patch.
data(0, dl + d2) <<
"," << ver[dl + d2][1]
4289 <<
">, <" << nrml[dl + d2][0] <<
", "
4290 << nrml[dl + d2][1] <<
", " << nrml[dl + d2][2]
4296 out <<
"triangle {" <<
'\n'
4297 <<
"\t<" << ver[dl][0] <<
"," << patch.
data(0, dl)
4298 <<
"," << ver[dl][1] <<
">," <<
'\n';
4299 out <<
"\t<" << ver[dl + d1][0] <<
","
4300 << patch.
data(0, dl + d1) <<
"," << ver[dl + d1][1]
4302 out <<
"\t<" << ver[dl + d1 + d2][0] <<
","
4303 << patch.
data(0, dl + d1 + d2) <<
","
4304 << ver[dl + d1 + d2][1] <<
">}" <<
'\n';
4307 out <<
"triangle {" <<
'\n'
4308 <<
"\t<" << ver[dl][0] <<
"," << patch.
data(0, dl)
4309 <<
"," << ver[dl][1] <<
">," <<
'\n';
4310 out <<
"\t<" << ver[dl + d1 + d2][0] <<
","
4311 << patch.
data(0, dl + d1 + d2) <<
","
4312 << ver[dl + d1 + d2][1] <<
">," <<
'\n';
4313 out <<
"\t<" << ver[dl + d2][0] <<
","
4314 << patch.
data(0, dl + d2) <<
"," << ver[dl + d2][1]
4322 Assert(n_subdivisions == 3,
4325 <<
"bicubic_patch {" <<
'\n'
4326 <<
" type 0" <<
'\n'
4327 <<
" flatness 0" <<
'\n'
4328 <<
" u_steps 0" <<
'\n'
4329 <<
" v_steps 0" <<
'\n';
4330 for (
int i = 0; i < 16; ++i)
4332 out <<
"\t<" << ver[i][0] <<
"," << patch.
data(0, i) <<
","
4333 << ver[i][1] <<
">";
4338 out <<
" texture {Tex}" <<
'\n' <<
"}" <<
'\n';
4342 if (!flags.bicubic_patch)
4345 out <<
" texture {Tex}" <<
'\n' <<
"}" <<
'\n' <<
'\n';
4357 template <
int dim,
int spacedim>
4361 const std::vector<std::string> &data_names,
4363 std::tuple<
unsigned int,
4370 do_write_povray(patches, data_names, flags, out);
4375 template <
int dim,
int spacedim>
4379 const std::vector<std::string> & ,
4381 std::tuple<
unsigned int,
4393 template <
int spacedim>
4397 const std::vector<std::string> & ,
4399 std::tuple<
unsigned int,
4408#ifndef DEAL_II_WITH_MPI
4417 if (patches.empty())
4427 std::multiset<EpsCell2d> cells;
4431 float min_color_value = std::numeric_limits<float>::max();
4432 float max_color_value = std::numeric_limits<float>::min();
4436 double heights[4] = {0, 0, 0, 0};
4440 for (
const auto &patch : patches)
4443 const unsigned int n = n_subdivisions + 1;
4444 const unsigned int d1 = 1;
4445 const unsigned int d2 = n;
4447 for (
unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
4448 for (
unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
4452 get_equispaced_location(patch, {i1, i2}, n_subdivisions);
4454 get_equispaced_location(patch, {i1 + 1, i2}, n_subdivisions);
4456 get_equispaced_location(patch, {i1, i2 + 1}, n_subdivisions);
4457 points[3] = get_equispaced_location(patch,
4465 patch.
data.n_rows() == 0,
4468 patch.
data.n_rows()));
4470 patch.
data.n_rows() != 0 ?
4474 heights[1] = patch.
data.n_rows() != 0 ?
4476 (i1 + 1) * d1 + i2 * d2) *
4479 heights[2] = patch.
data.n_rows() != 0 ?
4481 i1 * d1 + (i2 + 1) * d2) *
4484 heights[3] = patch.
data.n_rows() != 0 ?
4486 (i1 + 1) * d1 + (i2 + 1) * d2) *
4493 for (
unsigned int i = 0; i < 4; ++i)
4494 heights[i] = points[i][2];
4519 for (
unsigned int vertex = 0; vertex < 4; ++vertex)
4521 const double x = points[vertex][0], y = points[vertex][1],
4522 z = -heights[vertex];
4524 eps_cell.vertices[vertex][0] = -cz * x + sz * y;
4525 eps_cell.vertices[vertex][1] =
4526 -cx * sz * x - cx * cz * y - sx * z;
4550 (points[0] + points[1] + points[2] + points[3]) / 4;
4551 const double center_height =
4552 -(heights[0] + heights[1] + heights[2] + heights[3]) / 4;
4555 eps_cell.depth = -sx * sz * center_point[0] -
4556 sx * cz * center_point[1] + cx * center_height;
4561 patch.
data.n_rows() == 0,
4564 patch.
data.n_rows()));
4565 const double color_values[4] = {
4566 patch.
data.n_rows() != 0 ?
4570 patch.
data.n_rows() != 0 ?
4574 patch.
data.n_rows() != 0 ?
4578 patch.
data.n_rows() != 0 ?
4580 (i1 + 1) * d1 + (i2 + 1) * d2) :
4584 eps_cell.color_value = (color_values[0] + color_values[1] +
4585 color_values[3] + color_values[2]) /
4590 std::min(min_color_value, eps_cell.color_value);
4592 std::max(max_color_value, eps_cell.color_value);
4596 cells.insert(eps_cell);
4602 double x_min = cells.begin()->vertices[0][0];
4603 double x_max = x_min;
4604 double y_min = cells.begin()->vertices[0][1];
4605 double y_max = y_min;
4607 for (
const auto &cell : cells)
4608 for (
const auto &vertex : cell.vertices)
4610 x_min =
std::min(x_min, vertex[0]);
4611 x_max =
std::max(x_max, vertex[0]);
4612 y_min =
std::min(y_min, vertex[1]);
4613 y_max =
std::max(y_max, vertex[1]);
4618 const double scale =
4622 const Point<2> offset(x_min, y_min);
4627 out <<
"%!PS-Adobe-2.0 EPSF-1.2" <<
'\n'
4628 <<
"%%Title: deal.II Output" <<
'\n'
4629 <<
"%%Creator: the deal.II library" <<
'\n'
4632 <<
"%%BoundingBox: "
4636 <<
static_cast<unsigned int>((x_max - x_min) * scale + 0.5) <<
' '
4637 <<
static_cast<unsigned int>((y_max - y_min) * scale + 0.5) <<
'\n';
4646 out <<
"/m {moveto} bind def" <<
'\n'
4647 <<
"/l {lineto} bind def" <<
'\n'
4648 <<
"/s {setrgbcolor} bind def" <<
'\n'
4649 <<
"/sg {setgray} bind def" <<
'\n'
4650 <<
"/lx {lineto closepath stroke} bind def" <<
'\n'
4651 <<
"/lf {lineto closepath fill} bind def" <<
'\n';
4653 out <<
"%%EndProlog" <<
'\n' <<
'\n';
4655 out << flags.
line_width <<
" setlinewidth" <<
'\n';
4663 if (max_color_value == min_color_value)
4664 max_color_value = min_color_value + 1;
4668 for (
const auto &cell : cells)
4681 out << rgb_values.
red <<
" sg ";
4683 out << rgb_values.
red <<
' ' << rgb_values.
green <<
' '
4684 << rgb_values.
blue <<
" s ";
4689 out << (cell.vertices[0] - offset) * scale <<
" m "
4690 << (cell.vertices[1] - offset) * scale <<
" l "
4691 << (cell.vertices[3] - offset) * scale <<
" l "
4692 << (cell.vertices[2] - offset) * scale <<
" lf" <<
'\n';
4697 << (cell.vertices[0] - offset) * scale <<
" m "
4698 << (cell.vertices[1] - offset) * scale <<
" l "
4699 << (cell.vertices[3] - offset) * scale <<
" l "
4700 << (cell.vertices[2] - offset) * scale <<
" lx" <<
'\n';
4702 out <<
"showpage" <<
'\n';
4711 template <
int dim,
int spacedim>
4715 const std::vector<std::string> &data_names,
4717 std::tuple<
unsigned int,
4733#ifndef DEAL_II_WITH_MPI
4742 if (patches.empty())
4746 GmvStream gmv_out(out, flags);
4747 const unsigned int n_data_sets = data_names.size();
4750 Assert((patches[0].
data.n_rows() == n_data_sets &&
4751 !patches[0].points_are_available) ||
4752 (patches[0].data.n_rows() == n_data_sets + spacedim &&
4753 patches[0].points_are_available),
4755 (n_data_sets + spacedim) :
4757 patches[0].data.n_rows()));
4761 out <<
"gmvinput ascii" <<
'\n' <<
'\n';
4765 auto [n_nodes, n_cells] = count_nodes_and_cells(patches);
4780 [&patches]() {
return create_global_data_table(patches); });
4786 out <<
"nodes " << n_nodes <<
'\n';
4787 for (
unsigned int d = 0; d < spacedim; ++d)
4789 gmv_out.selected_component = d;
4795 for (
unsigned int d = spacedim; d < 3; ++d)
4797 for (
unsigned int i = 0; i < n_nodes; ++i)
4804 out <<
"cells " << n_cells <<
'\n';
4809 out <<
"variable" <<
'\n';
4813 std::move(*create_global_data_table_task.
return_value());
4817 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
4819 out << data_names[data_set] <<
" 1" <<
'\n';
4820 std::copy(data_vectors[data_set].
begin(),
4821 data_vectors[data_set].
end(),
4822 std::ostream_iterator<double>(out,
" "));
4823 out <<
'\n' <<
'\n';
4829 out <<
"endvars" <<
'\n';
4832 out <<
"endgmv" <<
'\n';
4843 template <
int dim,
int spacedim>
4847 const std::vector<std::string> &data_names,
4849 std::tuple<
unsigned int,
4863#ifndef DEAL_II_WITH_MPI
4872 if (patches.empty())
4876 TecplotStream tecplot_out(out, flags);
4878 const unsigned int n_data_sets = data_names.size();
4881 Assert((patches[0].
data.n_rows() == n_data_sets &&
4882 !patches[0].points_are_available) ||
4883 (patches[0].data.n_rows() == n_data_sets + spacedim &&
4884 patches[0].points_are_available),
4886 (n_data_sets + spacedim) :
4888 patches[0].data.n_rows()));
4892 auto [n_nodes, n_cells] = count_nodes_and_cells(patches);
4898 <<
"# This file was generated by the deal.II library." <<
'\n'
4902 <<
"# For a description of the Tecplot format see the Tecplot documentation."
4907 out <<
"Variables=";
4915 out <<
"\"x\", \"y\"";
4918 out <<
"\"x\", \"y\", \"z\"";
4924 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
4925 out <<
", \"" << data_names[data_set] <<
"\"";
4931 out <<
"t=\"" << flags.
zone_name <<
"\" ";
4934 out <<
"strandid=1, solutiontime=" << flags.
solution_time <<
", ";
4936 out <<
"f=feblock, n=" << n_nodes <<
", e=" << n_cells
4937 <<
", et=" << tecplot_cell_type[dim] <<
'\n';
4954 [&patches]() {
return create_global_data_table(patches); });
4960 for (
unsigned int d = 0; d < spacedim; ++d)
4962 tecplot_out.selected_component = d;
4973 std::move(*create_global_data_table_task.
return_value());
4976 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
4978 std::copy(data_vectors[data_set].
begin(),
4979 data_vectors[data_set].
end(),
4980 std::ostream_iterator<double>(out,
"\n"));
4995 template <
int dim,
int spacedim>
4999 const std::vector<std::string> &data_names,
5001 std::tuple<
unsigned int,
5005 &nonscalar_data_ranges,
5011#ifndef DEAL_II_WITH_MPI
5020 if (patches.empty())
5024 VtkStream vtk_out(out, flags);
5026 const unsigned int n_data_sets = data_names.size();
5028 if (patches[0].points_are_available)
5040 out <<
"# vtk DataFile Version 3.0" <<
'\n'
5041 <<
"#This file was generated by the deal.II library";
5049 out <<
'\n' <<
"ASCII" <<
'\n';
5051 out <<
"DATASET UNSTRUCTURED_GRID\n" <<
'\n';
5058 const unsigned int n_metadata =
5060 (flags.
time != std::numeric_limits<double>::lowest() ? 1 : 0));
5063 out <<
"FIELD FieldData " << n_metadata <<
'\n';
5067 out <<
"CYCLE 1 1 int\n" << flags.
cycle <<
'\n';
5069 if (flags.
time != std::numeric_limits<double>::lowest())
5071 out <<
"TIME 1 1 double\n" << flags.
time <<
'\n';
5077 unsigned int n_nodes;
5078 unsigned int n_cells;
5079 unsigned int n_points_and_n_cells;
5080 std::tie(n_nodes, n_cells, n_points_and_n_cells) =
5096 [&patches]() {
return create_global_data_table(patches); });
5102 out <<
"POINTS " << n_nodes <<
" double" <<
'\n';
5107 out <<
"CELLS " << n_cells <<
' ' << n_points_and_n_cells <<
'\n';
5115 out <<
"CELL_TYPES " << n_cells <<
'\n';
5119 for (
const auto &patch : patches)
5121 const auto vtk_cell_id =
5124 for (
unsigned int i = 0; i < vtk_cell_id[1]; ++i)
5125 out <<
' ' << vtk_cell_id[0];
5134 std::move(*create_global_data_table_task.
return_value());
5139 out <<
"POINT_DATA " << n_nodes <<
'\n';
5143 std::vector<bool> data_set_written(n_data_sets,
false);
5144 for (
const auto &nonscalar_data_range : nonscalar_data_ranges)
5149 "The VTK writer does not currently support outputting "
5150 "tensor data. Use the VTU writer instead."));
5153 std::get<0>(nonscalar_data_range),
5155 std::get<0>(nonscalar_data_range)));
5156 AssertThrow(std::get<1>(nonscalar_data_range) < n_data_sets,
5160 AssertThrow(std::get<1>(nonscalar_data_range) + 1 -
5161 std::get<0>(nonscalar_data_range) <=
5164 "Can't declare a vector with more than 3 components "
5168 for (
unsigned int i = std::get<0>(nonscalar_data_range);
5169 i <= std::get<1>(nonscalar_data_range);
5171 data_set_written[i] =
true;
5177 if (!std::get<2>(nonscalar_data_range).empty())
5178 out << std::get<2>(nonscalar_data_range);
5181 for (
unsigned int i = std::get<0>(nonscalar_data_range);
5182 i < std::get<1>(nonscalar_data_range);
5184 out << data_names[i] <<
"__";
5185 out << data_names[std::get<1>(nonscalar_data_range)];
5188 out <<
" double" <<
'\n';
5191 for (
unsigned int n = 0; n < n_nodes; ++n)
5193 switch (std::get<1>(nonscalar_data_range) -
5194 std::get<0>(nonscalar_data_range))
5197 out << data_vectors(std::get<0>(nonscalar_data_range), n)
5202 out << data_vectors(std::get<0>(nonscalar_data_range), n)
5204 << data_vectors(std::get<0>(nonscalar_data_range) + 1, n)
5208 out << data_vectors(std::get<0>(nonscalar_data_range), n)
5210 << data_vectors(std::get<0>(nonscalar_data_range) + 1, n)
5212 << data_vectors(std::get<0>(nonscalar_data_range) + 2, n)
5225 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
5226 if (data_set_written[data_set] ==
false)
5228 out <<
"SCALARS " << data_names[data_set] <<
" double 1" <<
'\n'
5229 <<
"LOOKUP_TABLE default" <<
'\n';
5230 std::copy(data_vectors[data_set].
begin(),
5231 data_vectors[data_set].
end(),
5232 std::ostream_iterator<double>(out,
" "));
5248 out <<
"<?xml version=\"1.0\" ?> \n";
5250 out <<
"# vtk DataFile Version 3.0" <<
'\n'
5251 <<
"#This file was generated by the deal.II library";
5262 out <<
"<VTKFile type=\"UnstructuredGrid\" version=\"2.2\"";
5264 out <<
"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\"";
5265 if (deal_ii_with_zlib &&
5267 out <<
" compressor=\"vtkZLibDataCompressor\"";
5268#ifdef DEAL_II_WORDS_BIGENDIAN
5269 out <<
" byte_order=\"BigEndian\"";
5271 out <<
" byte_order=\"LittleEndian\"";
5275 out <<
"<UnstructuredGrid>";
5285 out <<
" </UnstructuredGrid>\n";
5286 out <<
"</VTKFile>\n";
5291 template <
int dim,
int spacedim>
5295 const std::vector<std::string> &data_names,
5297 std::tuple<
unsigned int,
5301 &nonscalar_data_ranges,
5306 write_vtu_main(patches, data_names, nonscalar_data_ranges, flags, out);
5313 template <
int dim,
int spacedim>
5317 const std::vector<std::string> &data_names,
5319 std::tuple<
unsigned int,
5323 &nonscalar_data_ranges,
5337 unit.second.find(
'\"') == std::string::npos,
5339 "A physical unit you provided, <" + unit.second +
5340 ">, contained a quotation mark character. This is not allowed."));
5343#ifndef DEAL_II_WITH_MPI
5352 if (patches.empty())
5357 out <<
"<Piece NumberOfPoints=\"0\" NumberOfCells=\"0\" >\n"
5359 <<
"<DataArray type=\"UInt8\" Name=\"types\"></DataArray>\n"
5361 <<
" <PointData Scalars=\"scalars\">\n";
5362 std::vector<bool> data_set_written(data_names.size(),
false);
5363 for (
const auto &nonscalar_data_range : nonscalar_data_ranges)
5366 for (
unsigned int i = std::get<0>(nonscalar_data_range);
5367 i <= std::get<1>(nonscalar_data_range);
5369 data_set_written[i] =
true;
5373 out <<
" <DataArray type=\"Float32\" Name=\"";
5375 if (!std::get<2>(nonscalar_data_range).empty())
5376 out << std::get<2>(nonscalar_data_range);
5379 for (
unsigned int i = std::get<0>(nonscalar_data_range);
5380 i < std::get<1>(nonscalar_data_range);
5382 out << data_names[i] <<
"__";
5383 out << data_names[std::get<1>(nonscalar_data_range)];
5386 out <<
"\" NumberOfComponents=\"3\"></DataArray>\n";
5389 for (
unsigned int data_set = 0; data_set < data_names.size();
5391 if (data_set_written[data_set] ==
false)
5393 out <<
" <DataArray type=\"Float32\" Name=\""
5394 << data_names[data_set] <<
"\"></DataArray>\n";
5397 out <<
" </PointData>\n";
5398 out <<
"</Piece>\n";
5412 const unsigned int n_metadata =
5414 (flags.
time != std::numeric_limits<double>::lowest() ? 1 : 0));
5416 out <<
"<FieldData>\n";
5421 <<
"<DataArray type=\"Float32\" Name=\"CYCLE\" NumberOfTuples=\"1\" format=\"ascii\">"
5422 << flags.
cycle <<
"</DataArray>\n";
5424 if (flags.
time != std::numeric_limits<double>::lowest())
5427 <<
"<DataArray type=\"Float32\" Name=\"TIME\" NumberOfTuples=\"1\" format=\"ascii\">"
5428 << flags.
time <<
"</DataArray>\n";
5432 out <<
"</FieldData>\n";
5436 const unsigned int n_data_sets = data_names.size();
5439 if (patches[0].points_are_available)
5448 const char *ascii_or_binary =
5449 (deal_ii_with_zlib &&
5456 unsigned int n_nodes;
5457 unsigned int n_cells;
5458 std::tie(n_nodes, n_cells, std::ignore) =
5466 const auto stringize_vertex_information = [&patches,
5470 ascii_or_binary]() {
5471 std::ostringstream o;
5473 o <<
" <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\""
5474 << ascii_or_binary <<
"\">\n";
5475 const std::vector<Point<spacedim>> node_positions =
5480 std::vector<float> node_coordinates_3d;
5481 node_coordinates_3d.reserve(node_positions.size() * 3);
5482 for (
const auto &node_position : node_positions)
5484 for (
unsigned int d = 0; d < 3; ++d)
5486 node_coordinates_3d.emplace_back(node_position[d]);
5488 node_coordinates_3d.emplace_back(0.0f);
5490 vtu_stringize_array(node_coordinates_3d,
5495 o <<
" </DataArray>\n";
5496 o <<
" </Points>\n\n";
5505 const auto stringize_cell_to_vertex_information = [&patches,
5509 out.precision()]() {
5510 std::ostringstream o;
5513 o <<
" <DataArray type=\"Int32\" Name=\"connectivity\" format=\""
5514 << ascii_or_binary <<
"\">\n";
5516 std::vector<std::int32_t> cells;
5519 unsigned int first_vertex_of_patch = 0;
5521 std::vector<unsigned int> local_vertex_order;
5522 for (
const auto &patch : patches)
5524 local_vertex_order.clear();
5536 const unsigned int n_points = patch.
data.n_cols();
5537 Assert((dim == 2 && n_points == 6) ||
5538 (dim == 3 && n_points == 10),
5541 if (deal_ii_with_zlib &&
5545 for (
unsigned int i = 0; i < n_points; ++i)
5546 cells.push_back(first_vertex_of_patch + i);
5550 for (
unsigned int i = 0; i < n_points; ++i)
5551 o <<
'\t' << first_vertex_of_patch + i;
5555 first_vertex_of_patch += n_points;
5560 else if (patch.
reference_cell != ReferenceCells::get_hypercube<dim>())
5564 const unsigned int n_points = patch.
data.n_cols();
5566 if (deal_ii_with_zlib &&
5570 for (
unsigned int i = 0; i < n_points; ++i)
5572 first_vertex_of_patch +
5577 for (
unsigned int i = 0; i < n_points; ++i)
5579 << (first_vertex_of_patch +
5584 first_vertex_of_patch += n_points;
5589 const unsigned int n_points_per_direction = n_subdivisions + 1;
5592 const auto flush_current_cell = [&flags,
5595 first_vertex_of_patch,
5596 &local_vertex_order]() {
5597 if (deal_ii_with_zlib &&
5601 for (
const auto &c : local_vertex_order)
5602 cells.push_back(first_vertex_of_patch + c);
5606 for (
const auto &c : local_vertex_order)
5607 o <<
'\t' << first_vertex_of_patch + c;
5611 local_vertex_order.clear();
5616 local_vertex_order.reserve(Utilities::fixed_power<dim>(2));
5622 local_vertex_order.emplace_back(0);
5623 flush_current_cell();
5629 for (
unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
5631 const unsigned int starting_offset = i1;
5632 local_vertex_order.emplace_back(starting_offset);
5633 local_vertex_order.emplace_back(starting_offset +
5635 flush_current_cell();
5642 for (
unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
5643 for (
unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
5645 const unsigned int starting_offset =
5646 i2 * n_points_per_direction + i1;
5647 local_vertex_order.emplace_back(
5649 local_vertex_order.emplace_back(
5650 starting_offset + 1);
5651 local_vertex_order.emplace_back(
5652 starting_offset + n_points_per_direction + 1);
5653 local_vertex_order.emplace_back(
5654 starting_offset + n_points_per_direction);
5655 flush_current_cell();
5662 for (
unsigned int i3 = 0; i3 < n_subdivisions; ++i3)
5663 for (
unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
5664 for (
unsigned int i1 = 0; i1 < n_subdivisions;
5667 const unsigned int starting_offset =
5668 i3 * n_points_per_direction *
5669 n_points_per_direction +
5670 i2 * n_points_per_direction + i1;
5671 local_vertex_order.emplace_back(
5673 local_vertex_order.emplace_back(
5674 starting_offset + 1);
5675 local_vertex_order.emplace_back(
5676 starting_offset + n_points_per_direction +
5678 local_vertex_order.emplace_back(
5679 starting_offset + n_points_per_direction);
5680 local_vertex_order.emplace_back(
5681 starting_offset + n_points_per_direction *
5682 n_points_per_direction);
5683 local_vertex_order.emplace_back(
5685 n_points_per_direction *
5686 n_points_per_direction +
5688 local_vertex_order.emplace_back(
5690 n_points_per_direction *
5691 n_points_per_direction +
5692 n_points_per_direction + 1);
5693 local_vertex_order.emplace_back(
5695 n_points_per_direction *
5696 n_points_per_direction +
5697 n_points_per_direction);
5698 flush_current_cell();
5709 local_vertex_order.resize(
5710 Utilities::fixed_power<dim>(n_points_per_direction));
5712 if constexpr (dim == 0)
5716 "Point-like cells should not be possible "
5717 "when writing higher-order cells."));
5719 else if constexpr (dim == 1)
5721 for (
unsigned int i1 = 0; i1 < n_subdivisions + 1; ++i1)
5723 const unsigned int local_index = i1;
5724 const unsigned int connectivity_index =
5726 .vtk_lexicographic_to_node_index(
5730 local_vertex_order[connectivity_index] = local_index;
5732 flush_current_cell();
5734 else if constexpr (dim == 2)
5736 for (
unsigned int i2 = 0; i2 < n_subdivisions + 1; ++i2)
5737 for (
unsigned int i1 = 0; i1 < n_subdivisions + 1; ++i1)
5739 const unsigned int local_index =
5740 i2 * n_points_per_direction + i1;
5741 const unsigned int connectivity_index =
5743 .vtk_lexicographic_to_node_index(
5745 {{n_subdivisions, n_subdivisions}},
5747 local_vertex_order[connectivity_index] =
5750 flush_current_cell();
5752 else if constexpr (dim == 3)
5754 for (
unsigned int i3 = 0; i3 < n_subdivisions + 1; ++i3)
5755 for (
unsigned int i2 = 0; i2 < n_subdivisions + 1; ++i2)
5756 for (
unsigned int i1 = 0; i1 < n_subdivisions + 1;
5759 const unsigned int local_index =
5760 i3 * n_points_per_direction *
5761 n_points_per_direction +
5762 i2 * n_points_per_direction + i1;
5763 const unsigned int connectivity_index =
5765 .vtk_lexicographic_to_node_index(
5771 local_vertex_order[connectivity_index] =
5775 flush_current_cell();
5783 first_vertex_of_patch +=
5792 vtu_stringize_array(cells,
5798 o <<
" </DataArray>\n";
5818 const auto stringize_cell_offset_and_type_information =
5823 output_precision = out.precision()]() {
5824 std::ostringstream o;
5826 o <<
" <DataArray type=\"Int32\" Name=\"offsets\" format=\""
5827 << ascii_or_binary <<
"\">\n";
5829 std::vector<std::int32_t> offsets;
5830 offsets.reserve(n_cells);
5834 std::vector<unsigned int> cell_types;
5835 cell_types.reserve(n_cells);
5837 unsigned int first_vertex_of_patch = 0;
5839 for (
const auto &patch : patches)
5841 const auto vtk_cell_id =
5844 for (
unsigned int i = 0; i < vtk_cell_id[1]; ++i)
5846 cell_types.push_back(vtk_cell_id[0]);
5847 first_vertex_of_patch += vtk_cell_id[2];
5848 offsets.push_back(first_vertex_of_patch);
5852 vtu_stringize_array(offsets,
5857 o <<
" </DataArray>\n";
5859 o <<
" <DataArray type=\"UInt8\" Name=\"types\" format=\""
5860 << ascii_or_binary <<
"\">\n";
5862 if (deal_ii_with_zlib &&
5865 std::vector<std::uint8_t> cell_types_uint8_t(cell_types.size());
5866 for (
unsigned int i = 0; i < cell_types.size(); ++i)
5867 cell_types_uint8_t[i] =
static_cast<std::uint8_t
>(cell_types[i]);
5869 vtu_stringize_array(cell_types_uint8_t,
5876 vtu_stringize_array(cell_types,
5883 o <<
" </DataArray>\n";
5893 const auto stringize_nonscalar_data_range =
5899 output_precision = out.precision()](
const Table<2, float> &data_vectors,
5900 const auto &range) {
5901 std::ostringstream o;
5903 const auto first_component = std::get<0>(range);
5904 const auto last_component = std::get<1>(range);
5905 const auto &name = std::get<2>(range);
5906 const bool is_tensor =
5907 (std::get<3>(range) ==
5909 const unsigned int n_components = (is_tensor ? 9 : 3);
5916 AssertThrow((last_component + 1 - first_component <= 9),
5918 "Can't declare a tensor with more than 9 components "
5919 "in VTK/VTU format."));
5923 AssertThrow((last_component + 1 - first_component <= 3),
5925 "Can't declare a vector with more than 3 components "
5926 "in VTK/VTU format."));
5931 o <<
" <DataArray type=\"Float32\" Name=\"";
5937 for (
unsigned int i = first_component; i < last_component; ++i)
5938 o << data_names[i] <<
"__";
5939 o << data_names[last_component];
5942 o <<
"\" NumberOfComponents=\"" << n_components <<
"\" format=\""
5943 << ascii_or_binary <<
"\"";
5962 std::vector<float>
data;
5963 data.reserve(n_nodes * n_components);
5965 for (
unsigned int n = 0; n < n_nodes; ++n)
5969 switch (last_component - first_component)
5972 data.push_back(data_vectors(first_component, n));
5978 data.push_back(data_vectors(first_component, n));
5979 data.push_back(data_vectors(first_component + 1, n));
5984 data.push_back(data_vectors(first_component, n));
5985 data.push_back(data_vectors(first_component + 1, n));
5986 data.push_back(data_vectors(first_component + 2, n));
5999 const unsigned int size = last_component - first_component + 1;
6003 vtk_data[0][0] = data_vectors(first_component, n);
6008 for (
unsigned int c = 0; c <
size; ++c)
6012 vtk_data[ind[0]][ind[1]] =
6013 data_vectors(first_component + c, n);
6019 for (
unsigned int c = 0; c <
size; ++c)
6023 vtk_data[ind[0]][ind[1]] =
6024 data_vectors(first_component + c, n);
6035 for (
unsigned int i = 0; i < 3; ++i)
6036 for (
unsigned int j = 0; j < 3; ++j)
6037 data.push_back(vtk_data[i][j]);
6043 o <<
" </DataArray>\n";
6048 const auto stringize_scalar_data_set =
6052 output_precision = out.precision()](
const Table<2, float> &data_vectors,
6053 const unsigned int data_set) {
6054 std::ostringstream o;
6056 o <<
" <DataArray type=\"Float32\" Name=\"" << data_names[data_set]
6057 <<
"\" format=\"" << ascii_or_binary <<
"\"";
6061 o <<
" units=\"" << flags.
physical_units.at(data_names[data_set])
6066 const std::vector<float>
data(data_vectors[data_set].
begin(),
6067 data_vectors[data_set].
end());
6070 o <<
" </DataArray>\n";
6089 return create_global_data_table<dim, spacedim, float>(patches);
6103 std::move(*create_global_data_table_task.
return_value());
6109 std::vector<bool> data_set_handled(n_data_sets,
false);
6110 for (
const auto &range : nonscalar_data_ranges)
6113 const auto first_component = std::get<0>(range);
6114 const auto last_component = std::get<1>(range);
6115 for (
unsigned int i = first_component; i <= last_component; ++i)
6116 data_set_handled[i] =
true;
6119 return stringize_nonscalar_data_range(data_vectors, range);
6124 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
6125 if (data_set_handled[data_set] ==
false)
6128 return stringize_scalar_data_set(data_vectors, data_set);
6134 out <<
"<Piece NumberOfPoints=\"" << n_nodes <<
"\" NumberOfCells=\""
6135 << n_cells <<
"\" >\n";
6136 for (
const auto &s : mesh_tasks.return_values())
6138 out <<
" <PointData Scalars=\"scalars\">\n";
6141 out <<
" </PointData>\n";
6142 out <<
" </Piece>\n";
6156 const std::vector<std::string> &piece_names,
6157 const std::vector<std::string> &data_names,
6159 std::tuple<
unsigned int,
6163 &nonscalar_data_ranges,
6176 unit.second.find(
'\"') == std::string::npos,
6178 "A physical unit you provided, <" + unit.second +
6179 ">, contained a quotation mark character. This is not allowed."));
6182 const unsigned int n_data_sets = data_names.size();
6184 out <<
"<?xml version=\"1.0\"?>\n";
6187 out <<
"#This file was generated by the deal.II library"
6192 <<
"<VTKFile type=\"PUnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
6193 out <<
" <PUnstructuredGrid GhostLevel=\"0\">\n";
6201 const unsigned int n_metadata =
6203 (flags.
time != std::numeric_limits<double>::lowest() ? 1 : 0));
6205 out <<
" <FieldData>\n";
6210 <<
" <DataArray type=\"Float32\" Name=\"CYCLE\" NumberOfTuples=\"1\" format=\"ascii\">"
6211 << flags.
cycle <<
"</DataArray>\n";
6213 if (flags.
time != std::numeric_limits<double>::lowest())
6216 <<
" <DataArray type=\"Float32\" Name=\"TIME\" NumberOfTuples=\"1\" format=\"ascii\">"
6217 << flags.
time <<
"</DataArray>\n";
6221 out <<
" </FieldData>\n";
6224 out <<
" <PPointData Scalars=\"scalars\">\n";
6227 std::vector<bool> data_set_written(n_data_sets,
false);
6228 for (
const auto &nonscalar_data_range : nonscalar_data_ranges)
6230 const auto first_component = std::get<0>(nonscalar_data_range);
6231 const auto last_component = std::get<1>(nonscalar_data_range);
6232 const bool is_tensor =
6233 (std::get<3>(nonscalar_data_range) ==
6235 const unsigned int n_components = (is_tensor ? 9 : 3);
6242 AssertThrow((last_component + 1 - first_component <= 9),
6244 "Can't declare a tensor with more than 9 components "
6249 Assert((last_component + 1 - first_component <= 3),
6251 "Can't declare a vector with more than 3 components "
6256 for (
unsigned int i = std::get<0>(nonscalar_data_range);
6257 i <= std::get<1>(nonscalar_data_range);
6259 data_set_written[i] =
true;
6263 out <<
" <PDataArray type=\"Float32\" Name=\"";
6265 const std::string &name = std::get<2>(nonscalar_data_range);
6270 for (
unsigned int i = std::get<0>(nonscalar_data_range);
6271 i < std::get<1>(nonscalar_data_range);
6273 out << data_names[i] <<
"__";
6274 out << data_names[std::get<1>(nonscalar_data_range)];
6277 out <<
"\" NumberOfComponents=\"" << n_components
6278 <<
"\" format=\"ascii\"";
6290 data_names[std::get<1>(nonscalar_data_range)]) !=
6294 data_names[std::get<1>(nonscalar_data_range)])
6302 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
6303 if (data_set_written[data_set] ==
false)
6305 out <<
" <PDataArray type=\"Float32\" Name=\""
6306 << data_names[data_set] <<
"\" format=\"ascii\"";
6310 out <<
" units=\"" << flags.
physical_units.at(data_names[data_set])
6316 out <<
" </PPointData>\n";
6318 out <<
" <PPoints>\n";
6319 out <<
" <PDataArray type=\"Float32\" NumberOfComponents=\"3\"/>\n";
6320 out <<
" </PPoints>\n";
6322 for (
const auto &piece_name : piece_names)
6323 out <<
" <Piece Source=\"" << piece_name <<
"\"/>\n";
6325 out <<
" </PUnstructuredGrid>\n";
6326 out <<
"</VTKFile>\n";
6339 const std::vector<std::pair<double, std::string>> ×_and_names)
6343 out <<
"<?xml version=\"1.0\"?>\n";
6346 out <<
"#This file was generated by the deal.II library"
6351 <<
"<VTKFile type=\"Collection\" version=\"0.1\" ByteOrder=\"LittleEndian\">\n";
6352 out <<
" <Collection>\n";
6354 std::streamsize ss = out.precision();
6357 for (
const auto &time_and_name : times_and_names)
6358 out <<
" <DataSet timestep=\"" << time_and_name.first
6359 <<
"\" group=\"\" part=\"0\" file=\"" << time_and_name.second
6362 out <<
" </Collection>\n";
6363 out <<
"</VTKFile>\n";
6375 const std::vector<std::string> &piece_names)
6377 out <<
"!NBLOCKS " << piece_names.size() <<
'\n';
6378 for (
const auto &piece_name : piece_names)
6379 out << piece_name <<
'\n';
6388 const std::vector<std::vector<std::string>> &piece_names)
6392 if (piece_names.empty())
6395 const double nblocks = piece_names[0].size();
6397 ExcMessage(
"piece_names should be a vector of nonempty vectors."));
6399 out <<
"!NBLOCKS " << nblocks <<
'\n';
6400 for (
const auto &domain : piece_names)
6402 Assert(domain.size() == nblocks,
6404 "piece_names should be a vector of equal sized vectors."));
6405 for (
const auto &subdomain : domain)
6406 out << subdomain <<
'\n';
6417 const std::vector<std::pair<
double, std::vector<std::string>>>
6418 ×_and_piece_names)
6422 if (times_and_piece_names.empty())
6425 const double nblocks = times_and_piece_names[0].second.size();
6429 "time_and_piece_names should contain nonempty vectors of filenames for every timestep."));
6431 for (
const auto &domain : times_and_piece_names)
6432 out <<
"!TIME " << domain.first <<
'\n';
6434 out <<
"!NBLOCKS " << nblocks <<
'\n';
6435 for (
const auto &domain : times_and_piece_names)
6437 Assert(domain.second.size() == nblocks,
6439 "piece_names should be a vector of equal sized vectors."));
6440 for (
const auto &subdomain : domain.second)
6441 out << subdomain <<
'\n';
6449 template <
int dim,
int spacedim>
6453 const std::vector<std::string> &,
6455 std::tuple<
unsigned int,
6465 template <
int spacedim>
6469 const std::vector<std::string> & ,
6471 std::tuple<
unsigned int,
6480 unsigned int width = flags.
width;
6483 unsigned int margin_in_percent = 0;
6485 margin_in_percent = 5;
6489 double x_dimension, y_dimension, z_dimension;
6491 const auto &first_patch = patches[0];
6493 unsigned int n_subdivisions = first_patch.n_subdivisions;
6494 unsigned int n = n_subdivisions + 1;
6495 const unsigned int d1 = 1;
6496 const unsigned int d2 = n;
6499 std::array<Point<spacedim>, 4> projected_points;
6502 std::array<Point<2>, 4> projection_decompositions;
6505 get_equispaced_location(first_patch, {0, 0}, n_subdivisions);
6507 if (first_patch.data.n_rows() != 0)
6512 double x_min = projected_point[0];
6513 double x_max = x_min;
6514 double y_min = projected_point[1];
6515 double y_max = y_min;
6516 double z_min = first_patch.data.n_rows() != 0 ?
6519 double z_max = z_min;
6522 for (
const auto &patch : patches)
6525 n = n_subdivisions + 1;
6527 for (
unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
6529 for (
unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
6531 projected_points[0] =
6532 get_equispaced_location(patch, {i1, i2}, n_subdivisions);
6533 projected_points[1] =
6534 get_equispaced_location(patch, {i1 + 1, i2}, n_subdivisions);
6535 projected_points[2] =
6536 get_equispaced_location(patch, {i1, i2 + 1}, n_subdivisions);
6537 projected_points[3] = get_equispaced_location(patch,
6541 x_min =
std::min(x_min, projected_points[0][0]);
6542 x_min =
std::min(x_min, projected_points[1][0]);
6543 x_min =
std::min(x_min, projected_points[2][0]);
6544 x_min =
std::min(x_min, projected_points[3][0]);
6546 x_max =
std::max(x_max, projected_points[0][0]);
6547 x_max =
std::max(x_max, projected_points[1][0]);
6548 x_max =
std::max(x_max, projected_points[2][0]);
6549 x_max =
std::max(x_max, projected_points[3][0]);
6551 y_min =
std::min(y_min, projected_points[0][1]);
6552 y_min =
std::min(y_min, projected_points[1][1]);
6553 y_min =
std::min(y_min, projected_points[2][1]);
6554 y_min =
std::min(y_min, projected_points[3][1]);
6556 y_max =
std::max(y_max, projected_points[0][1]);
6557 y_max =
std::max(y_max, projected_points[1][1]);
6558 y_max =
std::max(y_max, projected_points[2][1]);
6559 y_max =
std::max(y_max, projected_points[3][1]);
6562 patch.
data.n_rows() == 0,
6565 patch.
data.n_rows()));
6567 z_min = std::min<double>(z_min,
6569 i1 * d1 + i2 * d2));
6570 z_min = std::min<double>(z_min,
6572 (i1 + 1) * d1 + i2 * d2));
6573 z_min = std::min<double>(z_min,
6575 i1 * d1 + (i2 + 1) * d2));
6577 std::min<double>(z_min,
6579 (i1 + 1) * d1 + (i2 + 1) * d2));
6581 z_max = std::max<double>(z_max,
6583 i1 * d1 + i2 * d2));
6584 z_max = std::max<double>(z_max,
6586 (i1 + 1) * d1 + i2 * d2));
6587 z_max = std::max<double>(z_max,
6589 i1 * d1 + (i2 + 1) * d2));
6591 std::max<double>(z_max,
6593 (i1 + 1) * d1 + (i2 + 1) * d2));
6598 x_dimension = x_max - x_min;
6599 y_dimension = y_max - y_min;
6600 z_dimension = z_max - z_min;
6607 float camera_focus = 0;
6610 camera_position[0] = 0.;
6611 camera_position[1] = 0.;
6612 camera_position[2] = z_min + 2. * z_dimension;
6614 camera_direction[0] = 0.;
6615 camera_direction[1] = 0.;
6616 camera_direction[2] = -1.;
6618 camera_horizontal[0] = 1.;
6619 camera_horizontal[1] = 0.;
6620 camera_horizontal[2] = 0.;
6622 camera_focus = .5 * z_dimension;
6628 const float angle_factor = 3.14159265f / 180.f;
6631 camera_position_temp[1] =
6634 camera_position_temp[2] =
6638 camera_direction_temp[1] =
6641 camera_direction_temp[2] =
6645 camera_horizontal_temp[1] =
6648 camera_horizontal_temp[2] =
6652 camera_position[1] = camera_position_temp[1];
6653 camera_position[2] = camera_position_temp[2];
6655 camera_direction[1] = camera_direction_temp[1];
6656 camera_direction[2] = camera_direction_temp[2];
6658 camera_horizontal[1] = camera_horizontal_temp[1];
6659 camera_horizontal[2] = camera_horizontal_temp[2];
6662 camera_position_temp[0] =
6665 camera_position_temp[1] =
6669 camera_direction_temp[0] =
6672 camera_direction_temp[1] =
6676 camera_horizontal_temp[0] =
6679 camera_horizontal_temp[1] =
6683 camera_position[0] = camera_position_temp[0];
6684 camera_position[1] = camera_position_temp[1];
6686 camera_direction[0] = camera_direction_temp[0];
6687 camera_direction[1] = camera_direction_temp[1];
6689 camera_horizontal[0] = camera_horizontal_temp[0];
6690 camera_horizontal[1] = camera_horizontal_temp[1];
6693 camera_position[0] = x_min + .5 * x_dimension;
6694 camera_position[1] = y_min + .5 * y_dimension;
6696 camera_position[0] += (z_min + 2. * z_dimension) *
6699 camera_position[1] -= (z_min + 2. * z_dimension) *
6705 double x_min_perspective, y_min_perspective;
6706 double x_max_perspective, y_max_perspective;
6707 double x_dimension_perspective, y_dimension_perspective;
6709 n_subdivisions = first_patch.n_subdivisions;
6710 n = n_subdivisions + 1;
6715 get_equispaced_location(first_patch, {0, 0}, n_subdivisions);
6717 if (first_patch.data.n_rows() != 0)
6722 point[0] = projected_point[0];
6723 point[1] = projected_point[1];
6724 point[2] = first_patch.data.n_rows() != 0 ?
6728 projection_decomposition = svg_project_point(point,
6734 x_min_perspective = projection_decomposition[0];
6735 x_max_perspective = projection_decomposition[0];
6736 y_min_perspective = projection_decomposition[1];
6737 y_max_perspective = projection_decomposition[1];
6740 for (
const auto &patch : patches)
6743 for (
unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
6745 for (
unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
6747 const std::array<Point<spacedim>, 4> projected_vertices{
6748 {get_equispaced_location(patch, {i1, i2}, n_subdivisions),
6749 get_equispaced_location(patch, {i1 + 1, i2}, n_subdivisions),
6750 get_equispaced_location(patch, {i1, i2 + 1}, n_subdivisions),
6751 get_equispaced_location(patch,
6756 patch.
data.n_rows() == 0,
6759 patch.
data.n_rows()));
6761 const std::array<Point<3>, 4> vertices = {
6762 {
Point<3>{projected_vertices[0][0],
6763 projected_vertices[0][1],
6764 patch.
data.n_rows() != 0 ?
6765 patch.
data(0, i1 * d1 + i2 * d2) :
6768 projected_vertices[1][1],
6769 patch.
data.n_rows() != 0 ?
6770 patch.
data(0, (i1 + 1) * d1 + i2 * d2) :
6773 projected_vertices[2][1],
6774 patch.
data.n_rows() != 0 ?
6775 patch.
data(0, i1 * d1 + (i2 + 1) * d2) :
6778 projected_vertices[3][1],
6779 patch.
data.n_rows() != 0 ?
6780 patch.
data(0, (i1 + 1) * d1 + (i2 + 1) * d2) :
6783 projection_decompositions = {
6784 {svg_project_point(vertices[0],
6789 svg_project_point(vertices[1],
6794 svg_project_point(vertices[2],
6799 svg_project_point(vertices[3],
6807 static_cast<double>(
6808 projection_decompositions[0][0]));
6811 static_cast<double>(
6812 projection_decompositions[1][0]));
6815 static_cast<double>(
6816 projection_decompositions[2][0]));
6819 static_cast<double>(
6820 projection_decompositions[3][0]));
6824 static_cast<double>(
6825 projection_decompositions[0][0]));
6828 static_cast<double>(
6829 projection_decompositions[1][0]));
6832 static_cast<double>(
6833 projection_decompositions[2][0]));
6836 static_cast<double>(
6837 projection_decompositions[3][0]));
6841 static_cast<double>(
6842 projection_decompositions[0][1]));
6845 static_cast<double>(
6846 projection_decompositions[1][1]));
6849 static_cast<double>(
6850 projection_decompositions[2][1]));
6853 static_cast<double>(
6854 projection_decompositions[3][1]));
6858 static_cast<double>(
6859 projection_decompositions[0][1]));
6862 static_cast<double>(
6863 projection_decompositions[1][1]));
6866 static_cast<double>(
6867 projection_decompositions[2][1]));
6870 static_cast<double>(
6871 projection_decompositions[3][1]));
6876 x_dimension_perspective = x_max_perspective - x_min_perspective;
6877 y_dimension_perspective = y_max_perspective - y_min_perspective;
6879 std::multiset<SvgCell> cells;
6882 for (
const auto &patch : patches)
6886 for (
unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
6888 for (
unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
6890 const std::array<Point<spacedim>, 4> projected_vertices = {
6891 {get_equispaced_location(patch, {i1, i2}, n_subdivisions),
6892 get_equispaced_location(patch, {i1 + 1, i2}, n_subdivisions),
6893 get_equispaced_location(patch, {i1, i2 + 1}, n_subdivisions),
6894 get_equispaced_location(patch,
6899 patch.
data.n_rows() == 0,
6902 patch.
data.n_rows()));
6906 cell.vertices[0][0] = projected_vertices[0][0];
6907 cell.vertices[0][1] = projected_vertices[0][1];
6908 cell.vertices[0][2] = patch.
data.n_rows() != 0 ?
6909 patch.
data(0, i1 * d1 + i2 * d2) :
6912 cell.vertices[1][0] = projected_vertices[1][0];
6913 cell.vertices[1][1] = projected_vertices[1][1];
6914 cell.vertices[1][2] = patch.
data.n_rows() != 0 ?
6915 patch.
data(0, (i1 + 1) * d1 + i2 * d2) :
6918 cell.vertices[2][0] = projected_vertices[2][0];
6919 cell.vertices[2][1] = projected_vertices[2][1];
6920 cell.vertices[2][2] = patch.
data.n_rows() != 0 ?
6921 patch.
data(0, i1 * d1 + (i2 + 1) * d2) :
6924 cell.vertices[3][0] = projected_vertices[3][0];
6925 cell.vertices[3][1] = projected_vertices[3][1];
6926 cell.vertices[3][2] =
6927 patch.
data.n_rows() != 0 ?
6928 patch.
data(0, (i1 + 1) * d1 + (i2 + 1) * d2) :
6931 cell.projected_vertices[0] =
6932 svg_project_point(cell.vertices[0],
6937 cell.projected_vertices[1] =
6938 svg_project_point(cell.vertices[1],
6943 cell.projected_vertices[2] =
6944 svg_project_point(cell.vertices[2],
6949 cell.projected_vertices[3] =
6950 svg_project_point(cell.vertices[3],
6956 cell.center = .25 * (cell.vertices[0] + cell.vertices[1] +
6957 cell.vertices[2] + cell.vertices[3]);
6958 cell.projected_center = svg_project_point(cell.center,
6964 cell.depth = cell.center.distance(camera_position);
6974 width =
static_cast<unsigned int>(
6975 .5 +
height * (x_dimension_perspective / y_dimension_perspective));
6976 unsigned int additional_width = 0;
6979 additional_width =
static_cast<unsigned int>(
6983 out <<
"<svg width=\"" << width + additional_width <<
"\" height=\""
6984 <<
height <<
"\" xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\">"
6986 <<
" <rect width=\"" << width + additional_width <<
"\" height=\""
6987 <<
height <<
"\" style=\"fill:white\"/>" <<
'\n'
6990 unsigned int triangle_counter = 0;
6993 for (
const auto &cell : cells)
6997 for (
unsigned int triangle_index = 0; triangle_index < 4;
7000 switch (triangle_index)
7003 points3d_triangle[0] = cell.vertices[0],
7004 points3d_triangle[1] = cell.vertices[1],
7005 points3d_triangle[2] = cell.center;
7008 points3d_triangle[0] = cell.vertices[1],
7009 points3d_triangle[1] = cell.vertices[3],
7010 points3d_triangle[2] = cell.center;
7013 points3d_triangle[0] = cell.vertices[3],
7014 points3d_triangle[1] = cell.vertices[2],
7015 points3d_triangle[2] = cell.center;
7018 points3d_triangle[0] = cell.vertices[2],
7019 points3d_triangle[1] = cell.vertices[0],
7020 points3d_triangle[2] = cell.center;
7027 svg_get_gradient_parameters(points3d_triangle);
7030 .667 - ((gradient_param[4] - z_min) / z_dimension) * .667;
7032 .667 - ((gradient_param[5] - z_min) / z_dimension) * .667;
7034 unsigned int start_r = 0;
7035 unsigned int start_g = 0;
7036 unsigned int start_b = 0;
7038 unsigned int stop_r = 0;
7039 unsigned int stop_g = 0;
7040 unsigned int stop_b = 0;
7042 unsigned int start_i =
static_cast<unsigned int>(start_h * 6.);
7043 unsigned int stop_i =
static_cast<unsigned int>(stop_h * 6.);
7045 double start_f = start_h * 6. - start_i;
7046 double start_q = 1. - start_f;
7048 double stop_f = stop_h * 6. - stop_i;
7049 double stop_q = 1. - stop_f;
7051 switch (start_i % 6)
7055 start_g =
static_cast<unsigned int>(.5 + 255. * start_f);
7058 start_r =
static_cast<unsigned int>(.5 + 255. * start_q),
7063 start_b =
static_cast<unsigned int>(.5 + 255. * start_f);
7066 start_g =
static_cast<unsigned int>(.5 + 255. * start_q),
7070 start_r =
static_cast<unsigned int>(.5 + 255. * start_f),
7075 start_b =
static_cast<unsigned int>(.5 + 255. * start_q);
7085 stop_g =
static_cast<unsigned int>(.5 + 255. * stop_f);
7088 stop_r =
static_cast<unsigned int>(.5 + 255. * stop_q),
7093 stop_b =
static_cast<unsigned int>(.5 + 255. * stop_f);
7096 stop_g =
static_cast<unsigned int>(.5 + 255. * stop_q),
7100 stop_r =
static_cast<unsigned int>(.5 + 255. * stop_f),
7105 stop_b =
static_cast<unsigned int>(.5 + 255. * stop_q);
7111 Point<3> gradient_start_point_3d, gradient_stop_point_3d;
7113 gradient_start_point_3d[0] = gradient_param[0];
7114 gradient_start_point_3d[1] = gradient_param[1];
7115 gradient_start_point_3d[2] = gradient_param[4];
7117 gradient_stop_point_3d[0] = gradient_param[2];
7118 gradient_stop_point_3d[1] = gradient_param[3];
7119 gradient_stop_point_3d[2] = gradient_param[5];
7122 svg_project_point(gradient_start_point_3d,
7128 svg_project_point(gradient_stop_point_3d,
7135 out <<
" <linearGradient id=\"" << triangle_counter
7136 <<
"\" gradientUnits=\"userSpaceOnUse\" "
7138 <<
static_cast<unsigned int>(
7140 ((gradient_start_point[0] - x_min_perspective) /
7141 x_dimension_perspective) *
7142 (width - (width / 100.) * 2. * margin_in_percent) +
7143 ((width / 100.) * margin_in_percent))
7146 <<
static_cast<unsigned int>(
7148 ((gradient_start_point[1] - y_min_perspective) /
7149 y_dimension_perspective) *
7153 <<
static_cast<unsigned int>(
7155 ((gradient_stop_point[0] - x_min_perspective) /
7156 x_dimension_perspective) *
7157 (width - (width / 100.) * 2. * margin_in_percent) +
7158 ((width / 100.) * margin_in_percent))
7161 <<
static_cast<unsigned int>(
7163 ((gradient_stop_point[1] - y_min_perspective) /
7164 y_dimension_perspective) *
7168 <<
" <stop offset=\"0\" style=\"stop-color:rgb(" << start_r
7169 <<
"," << start_g <<
"," << start_b <<
")\"/>" <<
'\n'
7170 <<
" <stop offset=\"1\" style=\"stop-color:rgb(" << stop_r
7171 <<
"," << stop_g <<
"," << stop_b <<
")\"/>" <<
'\n'
7172 <<
" </linearGradient>" <<
'\n';
7175 double x1 = 0, y1 = 0, x2 = 0, y2 = 0;
7176 double x3 = cell.projected_center[0];
7177 double y3 = cell.projected_center[1];
7179 switch (triangle_index)
7182 x1 = cell.projected_vertices[0][0],
7183 y1 = cell.projected_vertices[0][1],
7184 x2 = cell.projected_vertices[1][0],
7185 y2 = cell.projected_vertices[1][1];
7188 x1 = cell.projected_vertices[1][0],
7189 y1 = cell.projected_vertices[1][1],
7190 x2 = cell.projected_vertices[3][0],
7191 y2 = cell.projected_vertices[3][1];
7194 x1 = cell.projected_vertices[3][0],
7195 y1 = cell.projected_vertices[3][1],
7196 x2 = cell.projected_vertices[2][0],
7197 y2 = cell.projected_vertices[2][1];
7200 x1 = cell.projected_vertices[2][0],
7201 y1 = cell.projected_vertices[2][1],
7202 x2 = cell.projected_vertices[0][0],
7203 y2 = cell.projected_vertices[0][1];
7209 out <<
" <path d=\"M "
7210 <<
static_cast<unsigned int>(
7212 ((x1 - x_min_perspective) / x_dimension_perspective) *
7213 (width - (width / 100.) * 2. * margin_in_percent) +
7214 ((width / 100.) * margin_in_percent))
7216 <<
static_cast<unsigned int>(
7218 ((y1 - y_min_perspective) / y_dimension_perspective) *
7221 <<
static_cast<unsigned int>(
7223 ((x2 - x_min_perspective) / x_dimension_perspective) *
7224 (width - (width / 100.) * 2. * margin_in_percent) +
7225 ((width / 100.) * margin_in_percent))
7227 <<
static_cast<unsigned int>(
7229 ((y2 - y_min_perspective) / y_dimension_perspective) *
7232 <<
static_cast<unsigned int>(
7234 ((x3 - x_min_perspective) / x_dimension_perspective) *
7235 (width - (width / 100.) * 2. * margin_in_percent) +
7236 ((width / 100.) * margin_in_percent))
7238 <<
static_cast<unsigned int>(
7240 ((y3 - y_min_perspective) / y_dimension_perspective) *
7243 <<
static_cast<unsigned int>(
7245 ((x1 - x_min_perspective) / x_dimension_perspective) *
7246 (width - (width / 100.) * 2. * margin_in_percent) +
7247 ((width / 100.) * margin_in_percent))
7249 <<
static_cast<unsigned int>(
7251 ((y1 - y_min_perspective) / y_dimension_perspective) *
7253 <<
"\" style=\"stroke:black; fill:url(#" << triangle_counter
7264 out <<
'\n' <<
" <!-- colorbar -->" <<
'\n';
7266 unsigned int element_height =
static_cast<unsigned int>(
7267 ((
height / 100.) * (71. - 2. * margin_in_percent)) / 4);
7268 unsigned int element_width =
7269 static_cast<unsigned int>(.5 + (
height / 100.) * 2.5);
7271 additional_width = 0;
7274 static_cast<unsigned int>(.5 + (
height / 100.) * 2.5);
7276 for (
unsigned int index = 0; index < 4; ++index)
7278 double start_h = .667 - ((index + 1) / 4.) * .667;
7279 double stop_h = .667 - (index / 4.) * .667;
7281 unsigned int start_r = 0;
7282 unsigned int start_g = 0;
7283 unsigned int start_b = 0;
7285 unsigned int stop_r = 0;
7286 unsigned int stop_g = 0;
7287 unsigned int stop_b = 0;
7289 unsigned int start_i =
static_cast<unsigned int>(start_h * 6.);
7290 unsigned int stop_i =
static_cast<unsigned int>(stop_h * 6.);
7292 double start_f = start_h * 6. - start_i;
7293 double start_q = 1. - start_f;
7295 double stop_f = stop_h * 6. - stop_i;
7296 double stop_q = 1. - stop_f;
7298 switch (start_i % 6)
7302 start_g =
static_cast<unsigned int>(.5 + 255. * start_f);
7305 start_r =
static_cast<unsigned int>(.5 + 255. * start_q),
7310 start_b =
static_cast<unsigned int>(.5 + 255. * start_f);
7313 start_g =
static_cast<unsigned int>(.5 + 255. * start_q),
7317 start_r =
static_cast<unsigned int>(.5 + 255. * start_f),
7322 start_b =
static_cast<unsigned int>(.5 + 255. * start_q);
7332 stop_g =
static_cast<unsigned int>(.5 + 255. * stop_f);
7335 stop_r =
static_cast<unsigned int>(.5 + 255. * stop_q),
7340 stop_b =
static_cast<unsigned int>(.5 + 255. * stop_f);
7343 stop_g =
static_cast<unsigned int>(.5 + 255. * stop_q),
7347 stop_r =
static_cast<unsigned int>(.5 + 255. * stop_f),
7352 stop_b =
static_cast<unsigned int>(.5 + 255. * stop_q);
7359 out <<
" <linearGradient id=\"colorbar_" << index
7360 <<
"\" gradientUnits=\"userSpaceOnUse\" "
7361 <<
"x1=\"" << width + additional_width <<
"\" "
7363 <<
static_cast<unsigned int>(.5 + (
height / 100.) *
7364 (margin_in_percent + 29)) +
7365 (3 - index) * element_height
7367 <<
"x2=\"" << width + additional_width <<
"\" "
7369 <<
static_cast<unsigned int>(.5 + (
height / 100.) *
7370 (margin_in_percent + 29)) +
7371 (4 - index) * element_height
7374 <<
" <stop offset=\"0\" style=\"stop-color:rgb(" << start_r
7375 <<
"," << start_g <<
"," << start_b <<
")\"/>" <<
'\n'
7376 <<
" <stop offset=\"1\" style=\"stop-color:rgb(" << stop_r
7377 <<
"," << stop_g <<
"," << stop_b <<
")\"/>" <<
'\n'
7378 <<
" </linearGradient>" <<
'\n';
7383 <<
" x=\"" << width + additional_width <<
"\" y=\""
7384 <<
static_cast<unsigned int>(.5 + (
height / 100.) *
7385 (margin_in_percent + 29)) +
7386 (3 - index) * element_height
7387 <<
"\" width=\"" << element_width <<
"\" height=\""
7389 <<
"\" style=\"stroke:black; stroke-width:2; fill:url(#colorbar_"
7390 << index <<
")\"/>" <<
'\n';
7393 for (
unsigned int index = 0; index < 5; ++index)
7397 << width + additional_width +
7398 static_cast<unsigned int>(1.5 * element_width)
7400 <<
static_cast<unsigned int>(
7401 .5 + (
height / 100.) * (margin_in_percent + 29) +
7402 (4. - index) * element_height + 30.)
7404 <<
" style=\"text-anchor:start; font-size:80; font-family:Helvetica";
7406 if (index == 0 || index == 4)
7407 out <<
"; font-weight:bold";
7410 <<
static_cast<float>(
7411 (
static_cast<int>((z_min + index * (z_dimension / 4.)) *
7420 out <<
"</text>" <<
'\n';
7425 out <<
'\n' <<
"</svg>";
7431 template <
int dim,
int spacedim>
7435 const std::vector<std::string> &data_names,
7437 std::tuple<
unsigned int,
7441 &nonscalar_data_ranges,
7450 out << dim <<
' ' << spacedim <<
'\n';
7453 out <<
"[deal.II intermediate format graphics data]" <<
'\n'
7459 out << data_names.size() <<
'\n';
7460 for (
const auto &data_name : data_names)
7461 out << data_name <<
'\n';
7463 out << patches.size() <<
'\n';
7464 for (
unsigned int i = 0; i < patches.size(); ++i)
7465 out << patches[i] <<
'\n';
7467 out << nonscalar_data_ranges.size() <<
'\n';
7468 for (
const auto &nonscalar_data_range : nonscalar_data_ranges)
7469 out << std::get<0>(nonscalar_data_range) <<
' '
7470 << std::get<1>(nonscalar_data_range) <<
'\n'
7471 << std::get<2>(nonscalar_data_range) <<
'\n';
7479 template <
int dim,
int spacedim>
7483 const std::vector<std::string> &data_names,
7485 std::tuple<
unsigned int,
7489 &nonscalar_data_ranges,
7491 const std::string &filename,
7495#ifndef DEAL_II_WITH_MPI
7498 (void)nonscalar_data_ranges;
7505 ExcMessage(
"This functionality requires MPI to be enabled."));
7524 std::vector<char> my_buffer;
7526 boost::iostreams::filtering_ostream f;
7532# ifdef DEAL_II_WITH_ZLIB
7533 f.push(boost::iostreams::zlib_compressor(
7534 get_boost_zlib_compression_level(compression)));
7539 "Compression requires deal.II to be configured with ZLIB support."));
7542 boost::iostreams::back_insert_device<std::vector<char>> inserter(
7546 write_deal_II_intermediate<dim, spacedim>(
7547 patches, data_names, nonscalar_data_ranges, flags, f);
7549 const std::uint64_t my_size = my_buffer.size();
7555 const ParallelIntermediateHeader header{
7558 static_cast<std::uint64_t
>(compression),
7568 std::vector<std::uint64_t> chunk_sizes(n_ranks);
7569 int ierr = MPI_Gather(&my_size,
7571 Utilities::MPI::mpi_type_id_for_type<std::uint64_t>,
7572 static_cast<std::uint64_t *
>(chunk_sizes.data()),
7574 Utilities::MPI::mpi_type_id_for_type<std::uint64_t>,
7580 ierr = MPI_Info_create(&info);
7583 ierr = MPI_File_open(
7584 comm, filename.c_str(), MPI_MODE_CREATE | MPI_MODE_WRONLY, info, &fh);
7586 ierr = MPI_Info_free(&info);
7590 ierr = MPI_File_set_size(fh, 0);
7594 ierr = MPI_Barrier(
comm);
7601 fh, 0, &header,
sizeof(header), MPI_CHAR, MPI_STATUS_IGNORE);
7609 Utilities::MPI::mpi_type_id_for_type<std::uint64_t>,
7616 std::uint64_t prefix_sum = 0;
7617 ierr = MPI_Exscan(&my_size,
7620 Utilities::MPI::mpi_type_id_for_type<std::uint64_t>,
7626 const MPI_Offset offset =
static_cast<MPI_Offset
>(
sizeof(header)) +
7627 n_ranks *
sizeof(std::uint64_t) + prefix_sum;
7630 fh, offset, my_buffer.data(), my_size, MPI_CHAR, MPI_STATUS_IGNORE);
7639 ierr = MPI_File_sync(fh);
7642 ierr = MPI_File_close(&fh);
7649 std::pair<unsigned int, unsigned int>
7654 unsigned int dim, spacedim;
7655 input >> dim >> spacedim;
7657 return std::make_pair(dim, spacedim);
7666template <
int dim,
int spacedim>
7668 : default_subdivisions(1)
7674template <
int dim,
int spacedim>
7679 get_dataset_names(),
7680 get_nonscalar_data_ranges(),
7687template <
int dim,
int spacedim>
7692 get_dataset_names(),
7693 get_nonscalar_data_ranges(),
7700template <
int dim,
int spacedim>
7705 get_dataset_names(),
7706 get_nonscalar_data_ranges(),
7713template <
int dim,
int spacedim>
7718 get_dataset_names(),
7719 get_nonscalar_data_ranges(),
7726template <
int dim,
int spacedim>
7731 get_dataset_names(),
7732 get_nonscalar_data_ranges(),
7739template <
int dim,
int spacedim>
7744 get_dataset_names(),
7745 get_nonscalar_data_ranges(),
7752template <
int dim,
int spacedim>
7757 get_dataset_names(),
7758 get_nonscalar_data_ranges(),
7765template <
int dim,
int spacedim>
7770 get_dataset_names(),
7771 get_nonscalar_data_ranges(),
7776template <
int dim,
int spacedim>
7781 get_dataset_names(),
7782 get_nonscalar_data_ranges(),
7787template <
int dim,
int spacedim>
7792 get_dataset_names(),
7793 get_nonscalar_data_ranges(),
7799template <
int dim,
int spacedim>
7802 const std::string &filename,
7805#ifndef DEAL_II_WITH_MPI
7809 std::ofstream f(filename);
7817 int ierr = MPI_Info_create(&info);
7820 ierr = MPI_File_open(
7821 comm, filename.c_str(), MPI_MODE_CREATE | MPI_MODE_WRONLY, info, &fh);
7824 ierr = MPI_File_set_size(fh, 0);
7828 ierr = MPI_Barrier(
comm);
7830 ierr = MPI_Info_free(&info);
7834 unsigned int header_size;
7835 std::uint64_t footer_offset;
7840 std::stringstream ss;
7842 header_size = ss.str().size();
7845 fh, 0, ss.str().c_str(), header_size, MPI_CHAR, MPI_STATUS_IGNORE);
7849 ierr = MPI_Bcast(&header_size, 1, MPI_UNSIGNED, 0,
comm);
7853 const auto &patches = get_patches();
7862 std::stringstream ss;
7863 if (my_n_patches > 0 || (global_n_patches == 0 && myrank == 0))
7865 get_dataset_names(),
7866 get_nonscalar_data_ranges(),
7871 const std::uint64_t size_on_proc = ss.str().size();
7872 std::uint64_t prefix_sum = 0;
7873 ierr = MPI_Exscan(&size_on_proc,
7876 Utilities::MPI::mpi_type_id_for_type<std::uint64_t>,
7882 const MPI_Offset offset =
static_cast<MPI_Offset
>(header_size) + prefix_sum;
7892 if (myrank == n_ranks - 1)
7895 footer_offset = size_on_proc + offset;
7897 std::stringstream ss;
7899 const unsigned int footer_size = ss.str().size();
7917 ierr = MPI_File_sync(fh);
7920 ierr = MPI_File_close(&fh);
7927template <
int dim,
int spacedim>
7931 const std::vector<std::string> &piece_names)
const
7935 get_dataset_names(),
7936 get_nonscalar_data_ranges(),
7942template <
int dim,
int spacedim>
7945 const std::string &directory,
7946 const std::string &filename_without_extension,
7947 const unsigned int counter,
7949 const unsigned int n_digits_for_counter,
7950 const unsigned int n_groups)
const
7953 const unsigned int n_ranks =
7955 const unsigned int n_files_written =
7956 (n_groups == 0 || n_groups > n_ranks) ? n_ranks : n_groups;
7961 const unsigned int n_digits =
7964 const unsigned int color = rank % n_files_written;
7965 const std::string filename =
7966 directory + filename_without_extension +
"_" +
7970 if (n_groups == 0 || n_groups > n_ranks)
7973 std::ofstream output(filename);
7975 this->write_vtu(output);
7977 else if (n_groups == 1)
7980 this->write_vtu_in_parallel(filename, mpi_communicator);
7984#ifdef DEAL_II_WITH_MPI
7987 int ierr = MPI_Comm_split(mpi_communicator, color, rank, &comm_group);
7989 this->write_vtu_in_parallel(filename, comm_group);
7997 const std::string pvtu_filename =
7998 filename_without_extension +
"_" +
8003 std::vector<std::string> filename_vector;
8004 for (
unsigned int i = 0; i < n_files_written; ++i)
8006 const std::string filename =
8007 filename_without_extension +
"_" +
8011 filename_vector.emplace_back(filename);
8014 std::ofstream pvtu_output(directory + pvtu_filename);
8015 this->write_pvtu_record(pvtu_output, filename_vector);
8018 return pvtu_filename;
8023template <
int dim,
int spacedim>
8026 std::ostream &out)
const
8029 get_dataset_names(),
8030 get_nonscalar_data_ranges(),
8031 deal_II_intermediate_flags,
8037template <
int dim,
int spacedim>
8040 const std::string &filename,
8046 get_dataset_names(),
8047 get_nonscalar_data_ranges(),
8048 deal_II_intermediate_flags,
8056template <
int dim,
int spacedim>
8060 const std::string &h5_filename,
8061 const double cur_time,
8064 return create_xdmf_entry(
8065 data_filter, h5_filename, h5_filename, cur_time,
comm);
8070template <
int dim,
int spacedim>
8074 const std::string &h5_mesh_filename,
8075 const std::string &h5_solution_filename,
8076 const double cur_time,
8080 ExcMessage(
"XDMF only supports 2 or 3 space dimensions."));
8082#ifndef DEAL_II_WITH_HDF5
8086 (void)h5_mesh_filename;
8087 (void)h5_solution_filename;
8096 std::uint64_t local_node_cell_count[2], global_node_cell_count[2];
8098 local_node_cell_count[0] = data_filter.
n_nodes();
8099 local_node_cell_count[1] = data_filter.
n_cells();
8103 int ierr = MPI_Allreduce(local_node_cell_count,
8104 global_node_cell_count,
8106 Utilities::MPI::mpi_type_id_for_type<std::uint64_t>,
8118 const bool have_data = (data_filter.
n_nodes() > 0);
8121 const int key = myrank;
8122 const int color = (have_data ? 1 : 0);
8123 const int ierr = MPI_Comm_split(
comm, color, key, &split_comm);
8127 const bool am_i_first_rank_with_data =
8130 ierr = MPI_Comm_free(&split_comm);
8133 const int tag = 47381;
8136 if (am_i_first_rank_with_data)
8138 const auto &patches = get_patches();
8144 for (
const auto &patch : patches)
8150 h5_solution_filename,
8152 global_node_cell_count[0],
8153 global_node_cell_count[1],
8156 patches[0].reference_cell);
8157 const unsigned int n_data_sets = data_filter.
n_data_sets();
8161 for (
unsigned int i = 0; i < n_data_sets; ++i)
8171 ierr = MPI_Send(buffer.data(), buffer.size(), MPI_BYTE, 0, tag,
comm);
8180 if (myrank == 0 && !am_i_first_rank_with_data)
8185 int ierr = MPI_Probe(MPI_ANY_SOURCE, tag,
comm, &status);
8189 ierr = MPI_Get_count(&status, MPI_BYTE, &len);
8192 std::vector<char> buffer(len);
8193 ierr = MPI_Recv(buffer.data(),
8202 return Utilities::unpack<XDMFEntry>(buffer,
false);
8210template <
int dim,
int spacedim>
8213 const std::vector<XDMFEntry> &entries,
8214 const std::string &filename,
8217#ifdef DEAL_II_WITH_MPI
8221 const int myrank = 0;
8227 std::ofstream xdmf_file(filename);
8229 xdmf_file <<
"<?xml version=\"1.0\" ?>\n";
8230 xdmf_file <<
"<!DOCTYPE Xdmf SYSTEM \"Xdmf.dtd\" []>\n";
8231 xdmf_file <<
"<Xdmf Version=\"2.0\">\n";
8232 xdmf_file <<
" <Domain>\n";
8234 <<
" <Grid Name=\"CellTime\" GridType=\"Collection\" CollectionType=\"Temporal\">\n";
8236 for (
const auto &entry : entries)
8238 xdmf_file << entry.get_xdmf_content(3);
8241 xdmf_file <<
" </Grid>\n";
8242 xdmf_file <<
" </Domain>\n";
8243 xdmf_file <<
"</Xdmf>\n";
8255template <
int dim,
int spacedim>
8261 get_dataset_names(),
8262 get_nonscalar_data_ranges(),
8269#if defined(DEAL_II_WITH_HDF5) || defined(DEAL_II_WITH_NETCDF)
8273 struct DistributedMeshSizes
8275 uint64_t n_nodes_local;
8276 uint64_t n_cells_local;
8277 uint64_t n_nodes_global;
8278 uint64_t n_cells_global;
8279 uint64_t offset_nodes;
8280 uint64_t offset_cells;
8288 DistributedMeshSizes
8289 compute_global_mesh_size(
const uint64_t n_nodes_local,
8290 const uint64_t n_cells_local,
8293 const uint64_t n_local[2] = {n_nodes_local, n_cells_local};
8294 uint64_t n_global[2];
8295 uint64_t offsets[2] = {0, 0};
8297# ifdef DEAL_II_WITH_MPI
8299 MPI_Allreduce(n_local,
8315 n_global[0] = n_local[0];
8316 n_global[1] = n_local[1];
8322 n_local[0], n_local[1], n_global[0], n_global[1], offsets[0], offsets[1]};
8327#ifdef DEAL_II_WITH_HDF5
8331 template <
int dim,
int spacedim>
8336 const bool write_mesh_file,
8337 const std::string &mesh_filename,
8338 const std::string &solution_filename,
8341 hid_t h5_mesh_file_id = -1, h5_solution_file_id, file_plist_id, plist_id;
8342 hid_t node_dataspace, node_dataset, node_file_dataspace,
8343 node_memory_dataspace, node_dataset_id;
8344 hid_t cell_dataspace, cell_dataset, cell_file_dataspace,
8345 cell_memory_dataspace;
8346 hid_t pt_data_dataspace, pt_data_dataset, pt_data_file_dataspace,
8347 pt_data_memory_dataspace;
8349 hsize_t count[2], offset[2], node_ds_dim[2], cell_ds_dim[2];
8350 std::vector<double> node_data_vec;
8351 std::vector<unsigned int> cell_data_vec;
8354 file_plist_id = H5Pcreate(H5P_FILE_ACCESS);
8357# ifdef DEAL_II_WITH_MPI
8358# ifdef H5_HAVE_PARALLEL
8360 status = H5Pset_fapl_mpio(file_plist_id,
comm, MPI_INFO_NULL);
8365# ifndef DEAL_II_WITH_ZLIB
8371 const DistributedMeshSizes mesh_sizes =
8372 compute_global_mesh_size(data_filter.
n_nodes(),
8377 plist_id = H5Pcreate(H5P_DATASET_XFER);
8379# ifdef DEAL_II_WITH_MPI
8380# ifdef H5_HAVE_PARALLEL
8381 status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);
8386 if (write_mesh_file)
8389 h5_mesh_file_id = H5Fcreate(mesh_filename.c_str(),
8397 node_ds_dim[0] = mesh_sizes.n_nodes_global;
8398 node_ds_dim[1] = (spacedim < 2) ? 2 : spacedim;
8399 node_dataspace = H5Screate_simple(2, node_ds_dim,
nullptr);
8402 cell_ds_dim[0] = mesh_sizes.n_cells_global;
8403 cell_ds_dim[1] = patches[0].reference_cell.n_vertices();
8404 cell_dataspace = H5Screate_simple(2, cell_ds_dim,
nullptr);
8408# if H5Gcreate_vers == 1
8409 node_dataset = H5Dcreate(h5_mesh_file_id,
8415 node_dataset_id = H5Pcreate(H5P_DATASET_CREATE);
8416# ifdef DEAL_II_WITH_ZLIB
8417 H5Pset_deflate(node_dataset_id,
8419 H5Pset_chunk(node_dataset_id, 2, node_ds_dim);
8421 node_dataset = H5Dcreate(h5_mesh_file_id,
8428 H5Pclose(node_dataset_id);
8431# if H5Gcreate_vers == 1
8432 cell_dataset = H5Dcreate(h5_mesh_file_id,
8438 node_dataset_id = H5Pcreate(H5P_DATASET_CREATE);
8439# ifdef DEAL_II_WITH_ZLIB
8440 H5Pset_deflate(node_dataset_id,
8442 H5Pset_chunk(node_dataset_id, 2, cell_ds_dim);
8444 cell_dataset = H5Dcreate(h5_mesh_file_id,
8451 H5Pclose(node_dataset_id);
8456 status = H5Sclose(node_dataspace);
8458 status = H5Sclose(cell_dataspace);
8463 count[0] = mesh_sizes.n_nodes_local;
8464 count[1] = (spacedim < 2) ? 2 : spacedim;
8466 offset[0] = mesh_sizes.offset_nodes;
8469 node_memory_dataspace = H5Screate_simple(2, count,
nullptr);
8473 node_file_dataspace = H5Dget_space(node_dataset);
8475 status = H5Sselect_hyperslab(
8476 node_file_dataspace, H5S_SELECT_SET, offset,
nullptr, count,
nullptr);
8480 count[0] = mesh_sizes.n_cells_local;
8481 count[1] = patches[0].reference_cell.n_vertices();
8482 offset[0] = mesh_sizes.offset_cells;
8484 cell_memory_dataspace = H5Screate_simple(2, count,
nullptr);
8487 cell_file_dataspace = H5Dget_space(cell_dataset);
8489 status = H5Sselect_hyperslab(
8490 cell_file_dataspace, H5S_SELECT_SET, offset,
nullptr, count,
nullptr);
8495 status = H5Dwrite(node_dataset,
8497 node_memory_dataspace,
8498 node_file_dataspace,
8500 node_data_vec.data());
8502 node_data_vec.clear();
8505 data_filter.
fill_cell_data(mesh_sizes.offset_nodes, cell_data_vec);
8506 status = H5Dwrite(cell_dataset,
8508 cell_memory_dataspace,
8509 cell_file_dataspace,
8511 cell_data_vec.data());
8513 cell_data_vec.clear();
8516 status = H5Sclose(node_file_dataspace);
8518 status = H5Sclose(cell_file_dataspace);
8522 status = H5Sclose(node_memory_dataspace);
8524 status = H5Sclose(cell_memory_dataspace);
8528 status = H5Dclose(node_dataset);
8530 status = H5Dclose(cell_dataset);
8534 if (mesh_filename != solution_filename)
8536 status = H5Fclose(h5_mesh_file_id);
8542 if (mesh_filename == solution_filename && write_mesh_file)
8544 h5_solution_file_id = h5_mesh_file_id;
8549 h5_solution_file_id = H5Fcreate(solution_filename.c_str(),
8559 std::string vector_name;
8568 node_ds_dim[0] = mesh_sizes.n_nodes_global;
8569 node_ds_dim[1] = pt_data_vector_dim;
8570 pt_data_dataspace = H5Screate_simple(2, node_ds_dim,
nullptr);
8573# if H5Gcreate_vers == 1
8574 pt_data_dataset = H5Dcreate(h5_solution_file_id,
8575 vector_name.c_str(),
8580 node_dataset_id = H5Pcreate(H5P_DATASET_CREATE);
8581# ifdef DEAL_II_WITH_ZLIB
8582 H5Pset_deflate(node_dataset_id,
8584 H5Pset_chunk(node_dataset_id, 2, node_ds_dim);
8586 pt_data_dataset = H5Dcreate(h5_solution_file_id,
8587 vector_name.c_str(),
8593 H5Pclose(node_dataset_id);
8598 count[0] = mesh_sizes.n_nodes_local;
8599 count[1] = pt_data_vector_dim;
8600 offset[0] = mesh_sizes.offset_nodes;
8602 pt_data_memory_dataspace = H5Screate_simple(2, count,
nullptr);
8606 pt_data_file_dataspace = H5Dget_space(pt_data_dataset);
8608 status = H5Sselect_hyperslab(pt_data_file_dataspace,
8617 status = H5Dwrite(pt_data_dataset,
8619 pt_data_memory_dataspace,
8620 pt_data_file_dataspace,
8626 status = H5Sclose(pt_data_dataspace);
8628 status = H5Sclose(pt_data_memory_dataspace);
8630 status = H5Sclose(pt_data_file_dataspace);
8633 status = H5Dclose(pt_data_dataset);
8638 status = H5Pclose(file_plist_id);
8642 status = H5Pclose(plist_id);
8646 status = H5Fclose(h5_solution_file_id);
8654template <
int dim,
int spacedim>
8658 const std::vector<std::string> &data_names,
8660 std::tuple<
unsigned int,
8664 &nonscalar_data_ranges,
8667 const unsigned int n_data_sets = data_names.size();
8669#ifndef DEAL_II_WITH_MPI
8676 Assert(patches.size() > 0, ExcNoPatches());
8678 if (patches.empty())
8682 unsigned int n_nodes;
8683 std::tie(n_nodes, std::ignore) = count_nodes_and_cells(patches);
8698 [&patches]() {
return create_global_data_table(patches); });
8706 std::move(*create_global_data_table_task.
return_value());
8710 unsigned int i, n_th_vector, data_set, pt_data_vector_dim;
8711 std::string vector_name;
8712 for (n_th_vector = 0, data_set = 0; data_set < n_data_sets;)
8715 while (n_th_vector < nonscalar_data_ranges.size() &&
8716 std::get<0>(nonscalar_data_ranges[n_th_vector]) < data_set)
8720 if (n_th_vector < nonscalar_data_ranges.size() &&
8721 std::get<0>(nonscalar_data_ranges[n_th_vector]) == data_set)
8724 pt_data_vector_dim = std::get<1>(nonscalar_data_ranges[n_th_vector]) -
8725 std::get<0>(nonscalar_data_ranges[n_th_vector]) +
8730 std::get<1>(nonscalar_data_ranges[n_th_vector]) >=
8731 std::get<0>(nonscalar_data_ranges[n_th_vector]),
8732 ExcLowerRange(std::get<1>(nonscalar_data_ranges[n_th_vector]),
8733 std::get<0>(nonscalar_data_ranges[n_th_vector])));
8735 std::get<1>(nonscalar_data_ranges[n_th_vector]) < n_data_sets,
8736 ExcIndexRange(std::get<1>(nonscalar_data_ranges[n_th_vector]),
8742 if (!std::get<2>(nonscalar_data_ranges[n_th_vector]).empty())
8744 vector_name = std::get<2>(nonscalar_data_ranges[n_th_vector]);
8749 for (i = std::get<0>(nonscalar_data_ranges[n_th_vector]);
8750 i < std::get<1>(nonscalar_data_ranges[n_th_vector]);
8752 vector_name += data_names[i] +
"__";
8754 data_names[std::get<1>(nonscalar_data_ranges[n_th_vector])];
8760 pt_data_vector_dim = 1;
8761 vector_name = data_names[data_set];
8771 data_set += pt_data_vector_dim;
8777template <
int dim,
int spacedim>
8781 const std::string &filename,
8785 get_patches(), data_filter, hdf5_flags, filename,
comm);
8790template <
int dim,
int spacedim>
8794 const bool write_mesh_file,
8795 const std::string &mesh_filename,
8796 const std::string &solution_filename,
8810template <
int dim,
int spacedim>
8816 const std::string &filename,
8820 patches, data_filter, flags,
true, filename, filename,
comm);
8825template <
int dim,
int spacedim>
8831 const bool write_mesh_file,
8832 const std::string &mesh_filename,
8833 const std::string &solution_filename,
8839 "DataOutBase was asked to write HDF5 output for a space dimension of 1. "
8840 "HDF5 only supports datasets that live in 2 or 3 dimensions."));
8842#ifndef DEAL_II_WITH_HDF5
8848 (void)write_mesh_file;
8849 (void)mesh_filename;
8850 (void)solution_filename;
8859# ifndef H5_HAVE_PARALLEL
8863 "Serial HDF5 output on multiple processes is not yet supported."));
8873 Assert((patches.size() > 0) || (n_ranks > 1), ExcNoPatches());
8880 const bool have_patches = (patches.size() > 0);
8884 const int color = (have_patches ? 1 : 0);
8885 const int ierr = MPI_Comm_split(
comm, color, key, &split_comm);
8891 do_write_hdf5<dim, spacedim>(patches,
8900 const int ierr = MPI_Comm_free(&split_comm);
8908#ifdef DEAL_II_WITH_NETCDF
8921 create(
const std::string &filename,
8922 const int open_mode,
8926 auto ncerr = nc_create_par(
8927 filename.c_str(), open_mode,
comm, MPI_INFO_NULL, &ncfile.ncid);
8938 open(
const std::string &filename,
const int open_mode,
const MPI_Comm comm)
8941 auto ncerr = nc_open_par(
8942 filename.c_str(), open_mode,
comm, MPI_INFO_NULL, &ncfile.ncid);
8952 operator int()
const
8966 int ncerr = nc_close(ncid);
8968 ExcIO(
"Unexpected error when closing a NetCDF file."));
8991 nc_put_att_var(
const int ncid,
8997 [&](
const auto &v) {
8998 using T = std::remove_cv_t<std::remove_reference_t<
decltype(v)>>;
8999 if constexpr (std::is_same_v<T, std::string>)
9001 return nc_put_att_text(ncid, varid, name, v.size(), v.c_str());
9003 else if constexpr (std::is_same_v<T, double>)
9005 return nc_put_att_double(ncid, varid, name, NC_DOUBLE, 1, &v);
9007 else if constexpr (std::is_same_v<T, int>)
9009 return nc_put_att_int(ncid, varid, name, NC_INT, 1, &v);
9024 template <
class AttrMap = std::initializer_list<
9025 std::pair<std::string, DataOutBase::CFFlags::AttributeValue>>>
9027 write_cf_attributes(
int ncid,
int varid,
const AttrMap &attrs)
9029 for (
const auto &[att_name, att_value] : attrs)
9031 ncerr = nc_put_att_var(ncid, varid, att_name.c_str(), att_value);
9042 write_cf_user_defined_attributes(
int ncid,
9044 const std::string &var_name,
9047 if (
auto it_atts = flags.
attributes.find(var_name);
9050 write_cf_attributes(ncid, varid, it_atts->second);
9059 template <
int spacedim>
9062 std::vector<double> &x,
9063 std::vector<double> &y)
9065 std::vector<double> coords;
9067 x.resize(data_filter.
n_nodes());
9068 y.resize(data_filter.
n_nodes());
9069 for (
unsigned int i = 0; i < data_filter.
n_nodes(); ++i)
9071 x[i] = coords[i * spacedim];
9072 y[i] = coords[i * spacedim + 1];
9081 template <
int dim,
int spacedim>
9086 const std::string &filename,
9087 const DistributedMeshSizes &mesh_sizes,
9093 NcFile::create(filename, NC_WRITE | NC_CLOBBER | NC_NETCDF4,
comm);
9096 write_cf_attributes(
9100 {
"Conventions",
"CF-1.12"},
9102 "Created by the deal.II finite element library (https://dealii.org)."},
9104 write_cf_user_defined_attributes(ncid, flags,
"global", NC_GLOBAL);
9109 Assert(patches.size() > 0,
9113 Assert(std::all_of(patches.begin(),
9115 [ref = patches[0].reference_cell](
const auto &patch) {
9116 return patch.reference_cell == ref;
9119 "Mixed meshes are currently not supported in CF output."));
9123 int dim_node, dim_cell, dim_cell_node;
9124 ncerr = nc_def_dim(ncid,
"node", mesh_sizes.n_nodes_global, &dim_node);
9126 ncerr = nc_def_dim(ncid,
"face", mesh_sizes.n_cells_global, &dim_cell);
9128 int n_nodes_per_cell = patches[0].reference_cell.n_vertices();
9129 ncerr = nc_def_dim(ncid,
"face_node", n_nodes_per_cell, &dim_cell_node);
9135 int var_mesh, var_cell, var_x, var_y;
9136 ncerr = nc_def_var(ncid,
"mesh", NC_INT, 0,
nullptr, &var_mesh);
9138 write_cf_attributes(ncid,
9141 {
"cf_role",
"mesh_topology"},
9142 {
"long_name",
"Topology of a 2-d unstructured mesh"},
9143 {
"topology_dimension", 2},
9144 {
"node_coordinates",
"mesh_node_x mesh_node_y"},
9145 {
"face_node_connectivity",
"mesh_face_nodes"},
9146 {
"face_dimension",
"face"},
9150 int dummy_mesh_value = 0;
9151 ncerr = nc_put_var_int(ncid, var_mesh, &dummy_mesh_value);
9155 const int dims_cells[2] = {dim_cell, dim_cell_node};
9157 nc_def_var(ncid,
"mesh_face_nodes", NC_INT, 2, dims_cells, &var_cell);
9159 write_cf_attributes(ncid,
9162 {
"cf_role",
"face_node_connectivity"},
9163 {
"long_name",
"Corner nodes that make up each face."},
9171 ncerr = nc_def_var(ncid,
"mesh_node_x", NC_DOUBLE, 1, &dim_node, &var_x);
9173 write_cf_attributes(ncid, var_x, {{
"axis",
"X"}});
9174 write_cf_user_defined_attributes(ncid, flags,
"x", var_x);
9175 ncerr = nc_def_var(ncid,
"mesh_node_y", NC_DOUBLE, 1, &dim_node, &var_y);
9177 write_cf_attributes(ncid, var_y, {{
"axis",
"Y"}});
9178 write_cf_user_defined_attributes(ncid, flags,
"y", var_y);
9181 int dim_time, var_time;
9182 if (flags.
time > -std::numeric_limits<double>::infinity())
9184 ncerr = nc_def_dim(ncid,
"time", NC_UNLIMITED, &dim_time);
9186 ncerr = nc_def_var(ncid,
"time", NC_DOUBLE, 1, &dim_time, &var_time);
9188 write_cf_attributes(ncid,
9191 {
"standard_name",
"time"},
9193 {
"long_name",
"Time dimension"},
9195 write_cf_user_defined_attributes(ncid, flags,
"time", var_time);
9199 std::vector<double> x, y;
9200 get_node_coordinates<spacedim>(data_filter, x, y);
9201 const size_t coord_write_start[1] = {(size_t)mesh_sizes.offset_nodes};
9202 const size_t coord_write_count[1] = {(size_t)data_filter.
n_nodes()};
9203 ncerr = nc_put_vara_double(
9204 ncid, var_x, coord_write_start, coord_write_count, x.data());
9206 ncerr = nc_put_vara_double(
9207 ncid, var_y, coord_write_start, coord_write_count, y.data());
9211 std::vector<unsigned int> ucells;
9213 const size_t cell_write_start[2] = {(size_t)mesh_sizes.offset_cells,
9215 const size_t cell_write_count[2] = {(size_t)data_filter.
n_cells(),
9216 (size_t)n_nodes_per_cell};
9217 ncerr = nc_put_vara_uint(
9218 ncid, var_cell, cell_write_start, cell_write_count, ucells.data());
9228 write_cf_time_point(
int ncid,
int dim_time,
int var_time,
double time)
9236 ncerr = nc_var_par_access(ncid, var_time, NC_COLLECTIVE);
9240 ncerr = nc_inq_dimlen(ncid, dim_time, &n_time);
9242 auto times = std::vector<double>(n_time);
9243 ncerr = nc_get_var_double(ncid, var_time, times.data());
9245 const auto it_lb = std::lower_bound(times.begin(), times.end(), time);
9246 index_t = std::distance(times.begin(), it_lb);
9247 const size_t count_t = 1;
9248 ncerr = nc_put_vara_double(ncid, var_time, &index_t, &count_t, &time);
9251 ncerr = nc_var_par_access(ncid, var_time, NC_INDEPENDENT);
9265 const std::string &filename,
9266 const DistributedMeshSizes &mesh_sizes,
9272 const auto ncid = NcFile::open(filename, NC_WRITE | NC_NETCDF4,
comm);
9276 ncerr = nc_inq_dimid(ncid,
"node", &dim_node);
9280 int dim_time = -1, var_time = -1;
9282 const bool has_time = flags.
time > -std::numeric_limits<double>::infinity();
9285 int inq_dim_time = nc_inq_dimid(ncid,
"time", &dim_time);
9286 int inq_var_time = nc_inq_varid(ncid,
"time", &var_time);
9289 const bool has_time_dimension = inq_dim_time != NC_EBADDIM;
9290 const bool has_time_variable = inq_var_time != NC_ENOTVAR;
9291 AssertThrow(has_time_dimension && has_time_variable,
9293 "You are trying to add a time point to an existing file "
9294 "that does not have a time dimension or variable. "
9295 "Check the options in DataOutBase::cf_flags."));
9300 index_t = write_cf_time_point(ncid, dim_time, var_time, flags.
time);
9304 for (
auto dataset_index = 0u; dataset_index < data_filter.
n_data_sets();
9309 "Can only write scalar data in CF files."));
9316 const int inq_var = nc_inq_varid(ncid, key.c_str(), &varid);
9317 if (inq_var == NC_ENOTVAR)
9320 const int var_dims[2] = {dim_time, dim_node};
9321 ncerr = nc_def_var(ncid,
9325 &var_dims[has_time ? 0 : 1],
9328 write_cf_attributes(
9333 {
"location",
"node"},
9335 std::numeric_limits<double>::quiet_NaN()},
9338 write_cf_user_defined_attributes(ncid, flags, key, varid);
9347 ncerr = nc_var_par_access(ncid, varid, NC_COLLECTIVE);
9349 const size_t var_write_start[2] = {index_t,
9350 (size_t)mesh_sizes.offset_nodes};
9351 const size_t var_write_count[2] = {1, (size_t)data_filter.
n_nodes()};
9352 ncerr = nc_put_vara_double(
9355 &var_write_start[has_time ? 0 : 1],
9356 &var_write_count[has_time ? 0 : 1],
9359 ncerr = nc_var_par_access(ncid, varid, NC_INDEPENDENT);
9369template <
int dim,
int spacedim>
9373 const std::string &filename,
9377 get_patches(), data_filter, cf_flags, filename,
comm);
9382template <
int dim,
int spacedim>
9387 const std::string &filename,
9394 "DataOutBase can only write CF output in 2 dimensions."));
9397 "Mesh must have codimension 0 for writing CF output."));
9399#ifndef DEAL_II_WITH_NETCDF
9412 int ierr = MPI_Barrier(
comm);
9424 Assert((patches.size() > 0) || (n_ranks > 1), ExcNoPatches());
9431 const bool have_patches = (patches.size() > 0);
9435 const int color = (have_patches ? 1 : 0);
9436 const int ierr = MPI_Comm_split(
comm, color, key, &split_comm);
9439 ScopeExit cleanup_on_exit{[&split_comm] {
9440 int ierr = MPI_Comm_free(&split_comm);
9446 const DistributedMeshSizes mesh_sizes =
9447 compute_global_mesh_size(data_filter.
n_nodes(),
9455 mesh_sizes.n_nodes_global < uint64_t(std::numeric_limits<int>::max()),
9456 ExcMessage(
"The mesh has too many vertices to store in CF format."));
9462 patches, data_filter, flags, filename, mesh_sizes, split_comm);
9470 AssertThrow(std::filesystem::is_regular_file(filename),
9471 ExcIO(
"You are trying to write data to an existing file "
9472 "but no file of the given name exists. "
9473 "Create the file with a mesh by setting "
9474 "CFFlags::keep_existing_file = false."));
9479 write_cf_data(data_filter, flags, filename, mesh_sizes, split_comm);
9486template <
int dim,
int spacedim>
9494 output_format = default_fmt;
9496 switch (output_format)
9542 write_deal_II_intermediate(out);
9552template <
int dim,
int spacedim>
9561template <
int dim,
int spacedim>
9562template <
typename FlagType>
9566 if constexpr (std::is_same_v<FlagType, DataOutBase::DXFlags>)
9568 else if constexpr (std::is_same_v<FlagType, DataOutBase::UcdFlags>)
9570 else if constexpr (std::is_same_v<FlagType, DataOutBase::PovrayFlags>)
9571 povray_flags = flags;
9572 else if constexpr (std::is_same_v<FlagType, DataOutBase::EpsFlags>)
9574 else if constexpr (std::is_same_v<FlagType, DataOutBase::GmvFlags>)
9576 else if constexpr (std::is_same_v<FlagType, DataOutBase::Hdf5Flags>)
9578 else if constexpr (std::is_same_v<FlagType, DataOutBase::CFFlags>)
9580 else if constexpr (std::is_same_v<FlagType, DataOutBase::TecplotFlags>)
9581 tecplot_flags = flags;
9582 else if constexpr (std::is_same_v<FlagType, DataOutBase::VtkFlags>)
9584 else if constexpr (std::is_same_v<FlagType, DataOutBase::SvgFlags>)
9586 else if constexpr (std::is_same_v<FlagType, DataOutBase::GnuplotFlags>)
9587 gnuplot_flags = flags;
9588 else if constexpr (std::is_same_v<FlagType,
9590 deal_II_intermediate_flags = flags;
9597template <
int dim,
int spacedim>
9610template <
int dim,
int spacedim>
9617 "A name for the output format to be used");
9621 "Number of subdivisions of each mesh cell");
9667template <
int dim,
int spacedim>
9671 const std::string &output_name = prm.
get(
"Output format");
9673 default_subdivisions = prm.
get_integer(
"Subdivisions");
9676 dx_flags.parse_parameters(prm);
9680 ucd_flags.parse_parameters(prm);
9684 gnuplot_flags.parse_parameters(prm);
9688 povray_flags.parse_parameters(prm);
9692 eps_flags.parse_parameters(prm);
9696 gmv_flags.parse_parameters(prm);
9700 hdf5_flags.parse_parameters(prm);
9704 tecplot_flags.parse_parameters(prm);
9708 vtk_flags.parse_parameters(prm);
9712 deal_II_intermediate_flags.parse_parameters(prm);
9718template <
int dim,
int spacedim>
9722 return (
sizeof(default_fmt) +
9738template <
int dim,
int spacedim>
9740 std::tuple<
unsigned int,
9747 std::tuple<
unsigned int,
9754template <
int dim,
int spacedim>
9763 std::set<std::string> all_names;
9766 std::tuple<
unsigned int,
9770 ranges = this->get_nonscalar_data_ranges();
9771 const std::vector<std::string> data_names = this->get_dataset_names();
9772 const unsigned int n_data_sets = data_names.size();
9773 std::vector<bool> data_set_written(n_data_sets,
false);
9775 for (
const auto &range : ranges)
9777 const std::string &name = std::get<2>(range);
9780 Assert(all_names.find(name) == all_names.end(),
9782 "Error: names of fields in DataOut need to be unique, "
9784 name +
"' is used more than once."));
9785 all_names.insert(name);
9786 for (
unsigned int i = std::get<0>(range);
9787 i <= std::get<1>(range);
9789 data_set_written[i] =
true;
9793 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
9794 if (data_set_written[data_set] ==
false)
9796 const std::string &name = data_names[data_set];
9797 Assert(all_names.find(name) == all_names.end(),
9799 "Error: names of fields in DataOut need to be unique, "
9801 name +
"' is used more than once."));
9802 all_names.insert(name);
9812template <
int dim,
int spacedim>
9820 std::vector<typename ::DataOutBase::Patch<dim, spacedim>> tmp;
9824 std::vector<std::string> tmp;
9825 tmp.swap(dataset_names);
9829 std::tuple<
unsigned int,
9834 tmp.swap(nonscalar_data_ranges);
9841 std::pair<unsigned int, unsigned int> dimension_info =
9844 (dimension_info.second == spacedim),
9845 ExcIncompatibleDimensions(
9846 dimension_info.first, dim, dimension_info.second, spacedim));
9855 getline(in, header);
9857 std::ostringstream s;
9858 s <<
"[deal.II intermediate format graphics data]";
9860 Assert(header == s.str(), ExcUnexpectedInput(s.str(), header));
9864 getline(in, header);
9866 std::ostringstream s;
9870 Assert(header == s.str(), ExcUnexpectedInput(s.str(), header));
9874 getline(in, header);
9876 std::ostringstream s;
9880 Assert(header == s.str(),
9882 "Invalid or incompatible file format. Intermediate format "
9883 "files can only be read by the same deal.II version as they "
9884 "are written by."));
9888 unsigned int n_datasets;
9890 dataset_names.resize(n_datasets);
9891 for (
unsigned int i = 0; i < n_datasets; ++i)
9892 in >> dataset_names[i];
9894 unsigned int n_patches;
9896 patches.resize(n_patches);
9897 for (
unsigned int i = 0; i < n_patches; ++i)
9900 unsigned int n_nonscalar_data_ranges;
9901 in >> n_nonscalar_data_ranges;
9902 nonscalar_data_ranges.resize(n_nonscalar_data_ranges);
9903 for (
unsigned int i = 0; i < n_nonscalar_data_ranges; ++i)
9905 in >> std::get<0>(nonscalar_data_ranges[i]) >>
9906 std::get<1>(nonscalar_data_ranges[i]);
9915 std::get<2>(nonscalar_data_ranges[i]) = name;
9923template <
int dim,
int spacedim>
9929 ParallelIntermediateHeader header;
9930 in.read(
reinterpret_cast<char *
>(&header),
sizeof(header));
9932 header.magic == 0x00dea111,
9934 "Invalid header of parallel deal.II intermediate format encountered."));
9938 "Incorrect header version of parallel deal.II intermediate format."));
9940 std::vector<std::uint64_t> chunk_sizes(header.n_ranks);
9941 in.read(
reinterpret_cast<char *
>(chunk_sizes.data()),
9942 header.n_ranks *
sizeof(std::uint64_t));
9944 for (
unsigned int n = 0; n < header.n_ranks; ++n)
9948 std::vector<char> temp_buffer(chunk_sizes[n]);
9949 in.read(temp_buffer.data(), chunk_sizes[n]);
9952 header.compression) !=
9956 boost::iostreams::filtering_istreambuf f;
9959#ifdef DEAL_II_WITH_ZLIB
9960 f.push(boost::iostreams::zlib_decompressor());
9965 "Decompression requires deal.II to be configured with ZLIB support."));
9968 boost::iostreams::basic_array_source<char> source(temp_buffer.data(),
9969 temp_buffer.size());
9972 std::stringstream datastream;
9973 boost::iostreams::copy(f, datastream);
9983 temp_reader.
read(datastream);
9991template <
int dim,
int spacedim>
9995 using Patch = typename ::DataOutBase::Patch<dim, spacedim>;
9998 const std::vector<Patch> &source_patches = source.
get_patches();
10003 ExcIncompatibleDatasetNames());
10006 Assert(get_nonscalar_data_ranges().
size() ==
10008 ExcMessage(
"Both sources need to declare the same components "
10010 for (
unsigned int i = 0; i < get_nonscalar_data_ranges().size(); ++i)
10012 Assert(std::get<0>(get_nonscalar_data_ranges()[i]) ==
10014 ExcMessage(
"Both sources need to declare the same components "
10016 Assert(std::get<1>(get_nonscalar_data_ranges()[i]) ==
10018 ExcMessage(
"Both sources need to declare the same components "
10020 Assert(std::get<2>(get_nonscalar_data_ranges()[i]) ==
10022 ExcMessage(
"Both sources need to declare the same components "
10027 Assert(patches[0].n_subdivisions == source_patches[0].n_subdivisions,
10028 ExcIncompatiblePatchLists());
10029 Assert(patches[0].
data.n_rows() == source_patches[0].data.n_rows(),
10030 ExcIncompatiblePatchLists());
10031 Assert(patches[0].
data.n_cols() == source_patches[0].data.n_cols(),
10032 ExcIncompatiblePatchLists());
10036 const unsigned int old_n_patches = patches.size();
10037 patches.insert(patches.end(), source_patches.begin(), source_patches.end());
10040 for (
unsigned int i = old_n_patches; i < patches.size(); ++i)
10041 patches[i].patch_index += old_n_patches;
10044 for (
unsigned int i = old_n_patches; i < patches.size(); ++i)
10046 if (patches[i].neighbors[n] !=
10048 patches[i].neighbors[n] += old_n_patches;
10053template <
int dim,
int spacedim>
10054const std::vector<typename ::DataOutBase::Patch<dim, spacedim>> &
10062template <
int dim,
int spacedim>
10063std::vector<std::string>
10066 return dataset_names;
10071template <
int dim,
int spacedim>
10073 std::tuple<
unsigned int,
10079 return nonscalar_data_ranges;
10088 , h5_sol_filename(
"")
10089 , h5_mesh_filename(
"")
10091 , num_nodes(
numbers::invalid_unsigned_int)
10092 , num_cells(
numbers::invalid_unsigned_int)
10093 , dimension(
numbers::invalid_unsigned_int)
10094 , space_dimension(
numbers::invalid_unsigned_int)
10096 , n_vertices_per_cell(
numbers::invalid_unsigned_int)
10104 const std::uint64_t nodes,
10105 const std::uint64_t cells,
10106 const unsigned int dim_,
10108 :
XDMFEntry(filename, filename, time, nodes, cells, dim_, dim, cell_type)
10117 const std::string &solution_filename,
10119 const std::uint64_t nodes,
10120 const std::uint64_t cells,
10121 const unsigned int dim_,
10147 if (cell_type == ReferenceCells::Invalid<dim>)
10148 return ReferenceCells::get_hypercube<dim>();
10158 const std::string &solution_filename,
10160 const std::uint64_t nodes,
10161 const std::uint64_t cells,
10162 const unsigned int dim_,
10163 const unsigned int spacedim,
10166 , h5_sol_filename(solution_filename)
10167 , h5_mesh_filename(mesh_filename)
10172 , space_dimension(spacedim)
10177 cell_type_hex_if_invalid<dim>(cell_type_);
10180 if constexpr (dim == 0)
10184 else if constexpr (dim == 1)
10188 else if constexpr (dim == 2)
10203 else if constexpr (dim == 3)
10226 const unsigned int dimension)
10239 indent(
const unsigned int indent_level)
10241 std::string res =
"";
10242 for (
unsigned int i = 0; i < indent_level; ++i)
10256 std::stringstream ss;
10258 ss << indent(indent_level + 0)
10259 <<
"<Grid Name=\"mesh\" GridType=\"Uniform\">\n";
10260 ss << indent(indent_level + 1) <<
"<Time Value=\"" <<
entry_time <<
"\"/>\n";
10261 ss << indent(indent_level + 1) <<
"<Geometry GeometryType=\""
10263 ss << indent(indent_level + 2) <<
"<DataItem Dimensions=\"" <<
num_nodes
10265 <<
"\" NumberType=\"Float\" Precision=\"8\" Format=\"HDF\">\n";
10267 ss << indent(indent_level + 2) <<
"</DataItem>\n";
10268 ss << indent(indent_level + 1) <<
"</Geometry>\n";
10273 ss << indent(indent_level + 1) <<
"<Topology TopologyType=\""
10276 ss <<
"\" NodesPerElement=\"1\">\n";
10278 ss <<
"\" NodesPerElement=\"2\">\n";
10283 ss << indent(indent_level + 2) <<
"<DataItem Dimensions=\"" <<
num_cells
10285 <<
"\" NumberType=\"UInt\" Format=\"HDF\">\n";
10288 ss << indent(indent_level + 2) <<
"</DataItem>\n";
10289 ss << indent(indent_level + 1) <<
"</Topology>\n";
10295 ss << indent(indent_level + 1)
10296 <<
"<Topology TopologyType=\"Polyvertex\" NumberOfElements=\""
10298 ss << indent(indent_level + 1) <<
"</Topology>\n";
10303 ss << indent(indent_level + 1) <<
"<Attribute Name=\""
10304 << attribute_dim.first <<
"\" AttributeType=\""
10305 << (attribute_dim.second > 1 ?
"Vector" :
"Scalar")
10306 <<
"\" Center=\"Node\">\n";
10308 ss << indent(indent_level + 2) <<
"<DataItem Dimensions=\"" <<
num_nodes
10309 <<
" " << (attribute_dim.second > 1 ? 3 : 1)
10310 <<
"\" NumberType=\"Float\" Precision=\"8\" Format=\"HDF\">\n";
10312 << attribute_dim.first <<
'\n';
10313 ss << indent(indent_level + 2) <<
"</DataItem>\n";
10314 ss << indent(indent_level + 1) <<
"</Attribute>\n";
10317 ss << indent(indent_level + 0) <<
"</Grid>\n";
10326 template <
int dim,
int spacedim>
10331 out <<
"[deal.II intermediate Patch<" << dim <<
',' << spacedim <<
">]"
10338 for (
const unsigned int i : patch.
reference_cell.vertex_indices())
10342 for (
const unsigned int i : patch.
reference_cell.face_indices())
10350 out << patch.
data.n_rows() <<
' ' << patch.
data.n_cols() <<
'\n';
10351 for (
unsigned int i = 0; i < patch.
data.n_rows(); ++i)
10352 for (
unsigned int j = 0; j < patch.
data.n_cols(); ++j)
10353 out << patch.
data[i][j] <<
' ';
10362 template <
int dim,
int spacedim>
10371 std::string header;
10374 getline(in, header);
10375 while ((header.size() != 0) && (header.back() ==
' '))
10376 header.erase(header.size() - 1);
10378 while ((header.empty()) && in);
10380 std::ostringstream s;
10381 s <<
"[deal.II intermediate Patch<" << dim <<
',' << spacedim <<
">]";
10383 Assert(header == s.str(), ExcUnexpectedInput(s.str(), header));
10387 if constexpr (dim > 0)
10391 for (
const unsigned int i : patch.
reference_cell.vertex_indices())
10394 for (
const unsigned int i : patch.
reference_cell.face_indices())
10401 unsigned int n_subdivisions;
10402 in >> n_subdivisions;
10403 if constexpr (dim > 1)
10408 unsigned int n_rows, n_cols;
10409 in >> n_rows >> n_cols;
10410 patch.
data.reinit(n_rows, n_cols);
10411 for (
unsigned int i = 0; i < patch.
data.n_rows(); ++i)
10412 for (
unsigned int j = 0; j < patch.
data.n_cols(); ++j)
10413 in >> patch.
data[i][j];
10424#include "base/data_out_base.inst"
const double * get_data_set(const unsigned int set_num) const
std::map< unsigned int, unsigned int > filtered_points
std::string get_data_set_name(const unsigned int set_num) const
void internal_add_cell(const unsigned int cell_index, const unsigned int pt_index)
void write_cell(const unsigned int index, const unsigned int start, const std::array< unsigned int, dim > &offsets)
std::vector< unsigned int > data_set_dims
unsigned int n_nodes() const
void fill_node_data(std::vector< double > &node_data) const
std::vector< std::vector< double > > data_sets
unsigned int get_data_set_dim(const unsigned int set_num) const
void write_data_set(const std::string &name, const unsigned int dimension, const unsigned int set_num, const Table< 2, double > &data_vectors)
std::map< unsigned int, unsigned int > filtered_cells
std::vector< std::string > data_set_names
void fill_cell_data(const unsigned int local_node_offset, std::vector< unsigned int > &cell_data) const
void write_point(const unsigned int index, const Point< dim > &p)
void write_cell_single(const unsigned int index, const unsigned int start, const unsigned int n_points, const ReferenceCell< dim > &reference_cell)
unsigned int n_cells() const
DataOutBase::DataOutFilterFlags flags
Map3DPoint existing_points
unsigned int n_data_sets() const
XDMFEntry create_xdmf_entry(const DataOutBase::DataOutFilter &data_filter, const std::string &h5_filename, const double cur_time, const MPI_Comm comm) const
virtual std::vector< std::tuple< unsigned int, unsigned int, std::string, DataComponentInterpretation::DataComponentInterpretation > > get_nonscalar_data_ranges() const
void parse_parameters(ParameterHandler &prm)
void write_filtered_data(DataOutBase::DataOutFilter &filtered_data) const
void write_cf_parallel(const DataOutBase::DataOutFilter &data_filter, const std::string &filename, const MPI_Comm comm) const
void write_pvtu_record(std::ostream &out, const std::vector< std::string > &piece_names) const
static void declare_parameters(ParameterHandler &prm)
void write_ucd(std::ostream &out) const
void write_povray(std::ostream &out) const
std::string default_suffix(const DataOutBase::OutputFormat output_format=DataOutBase::default_format) const
void write_xdmf_file(const std::vector< XDMFEntry > &entries, const std::string &filename, const MPI_Comm comm) const
std::size_t memory_consumption() const
void set_default_format(const DataOutBase::OutputFormat default_format)
void write(std::ostream &out, const DataOutBase::OutputFormat output_format=DataOutBase::default_format) const
void write_gnuplot(std::ostream &out) const
void write_vtu(std::ostream &out) const
void write_hdf5_parallel(const DataOutBase::DataOutFilter &data_filter, const std::string &filename, const MPI_Comm comm) const
void write_tecplot(std::ostream &out) const
void write_deal_II_intermediate_in_parallel(const std::string &filename, const MPI_Comm comm, const DataOutBase::CompressionLevel compression) const
void write_svg(std::ostream &out) const
void write_vtu_in_parallel(const std::string &filename, const MPI_Comm comm) const
void validate_dataset_names() const
void set_flags(const FlagType &flags)
void write_vtk(std::ostream &out) const
void write_gmv(std::ostream &out) const
std::string write_vtu_with_pvtu_record(const std::string &directory, const std::string &filename_without_extension, const unsigned int counter, const MPI_Comm mpi_communicator, const unsigned int n_digits_for_counter=numbers::invalid_unsigned_int, const unsigned int n_groups=0) const
void write_eps(std::ostream &out) const
void write_dx(std::ostream &out) const
void write_deal_II_intermediate(std::ostream &out) const
void merge(const DataOutReader< dim, spacedim > &other)
void read_whole_parallel_file(std::istream &in)
virtual const std::vector<::DataOutBase::Patch< dim, spacedim > > & get_patches() const override
void read(std::istream &in)
virtual std::vector< std::tuple< unsigned int, unsigned int, std::string, DataComponentInterpretation::DataComponentInterpretation > > get_nonscalar_data_ranges() const override
virtual std::vector< std::string > get_dataset_names() const override
void enter_subsection(const std::string &subsection, const bool create_path_if_needed=true)
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
constexpr unsigned int n_vertices() const
static constexpr TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
std::vector< RT > return_values()
internal::return_value< RT >::reference_type return_value()
unsigned int n_vertices_per_cell
std::string h5_sol_filename
std::string cell_type_name
std::string h5_mesh_filename
std::string get_xdmf_content(const unsigned int indent_level) const
void add_attribute(const std::string &attr_name, const unsigned int dimension)
std::map< std::string, unsigned int > attribute_dims
unsigned int space_dimension
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_PACKAGE_VERSION
constexpr bool running_in_debug_mode()
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_PACKAGE_NAME
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
#define AssertThrowNC(code)
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcFileNotOpen(std::string arg1)
static ::ExceptionBase & ExcNotEnoughSpaceDimensionLabels()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcNoPatches()
static ::ExceptionBase & ExcNeedsNetCDF()
#define DeclException2(Exception2, type1, type2, outsequence)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcLowerRange(int arg1, int arg2)
#define AssertThrowMPI(error_code)
#define AssertNothrow(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNeedsHDF5()
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcInvalidDatasetSize(int arg1, int arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Task< RT > new_task(const std::function< RT()> &function)
const unsigned int my_rank
std::vector< index_type > data
DataComponentInterpretation
@ component_is_part_of_tensor
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)
std::pair< unsigned int, unsigned int > determine_intermediate_format_dimensions(std::istream &input)
std::ostream & operator<<(std::ostream &out, const Patch< dim, spacedim > &patch)
void write_nodes(const std::vector< Patch< dim, spacedim > > &patches, StreamType &out)
void write_deal_II_intermediate_in_parallel(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 Deal_II_IntermediateFlags &flags, const std::string &filename, const MPI_Comm comm, const CompressionLevel compression)
void write_ucd(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 UcdFlags &flags, std::ostream &out)
void write_cf_parallel(const std::vector< Patch< dim, spacedim > > &patches, const DataOutFilter &data_filter, const DataOutBase::CFFlags &flags, const std::string &filename, const MPI_Comm comm)
void write_dx(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 DXFlags &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_gmv(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 GmvFlags &flags, std::ostream &out)
void write_data(const std::vector< Patch< dim, spacedim > > &patches, unsigned int n_data_sets, const bool double_precision, StreamType &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_pvd_record(std::ostream &out, const std::vector< std::pair< double, std::string > > ×_and_names)
void write_deal_II_intermediate(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 Deal_II_IntermediateFlags &flags, std::ostream &out)
void write_vtu_footer(std::ostream &out)
void write_cells(const std::vector< Patch< dim, spacedim > > &patches, StreamType &out)
void write_tecplot(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 TecplotFlags &flags, std::ostream &out)
void write_filtered_data(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, DataOutFilter &filtered_data)
OutputFormat parse_output_format(const std::string &format_name)
std::vector< Point< spacedim > > get_node_positions(const std::vector< Patch< dim, spacedim > > &patches)
void write_hdf5_parallel(const std::vector< Patch< dim, spacedim > > &patches, const DataOutFilter &data_filter, const DataOutBase::Hdf5Flags &flags, const std::string &filename, const MPI_Comm comm)
std::istream & operator>>(std::istream &in, Patch< dim, spacedim > &patch)
std::string get_output_format_names()
void write_svg(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 SvgFlags &flags, std::ostream &out)
void write_visit_record(std::ostream &out, const std::vector< std::string > &piece_names)
void write_high_order_cells(const std::vector< Patch< dim, spacedim > > &patches, StreamType &out, const bool legacy_format)
std::string default_suffix(const OutputFormat output_format)
void write_povray(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 PovrayFlags &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_vtk(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_pvtu_record(std::ostream &out, const std::vector< std::string > &piece_names, 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)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr ReferenceCell< 3 > Hexahedron
constexpr ReferenceCell< 2 > Quadrilateral
constexpr ReferenceCell< 1 > Line
constexpr ReferenceCell< 2 > Triangle
constexpr ReferenceCell< 3 > Tetrahedron
constexpr ReferenceCell< 3 > Pyramid
constexpr ReferenceCell< 3 > Wedge
constexpr ReferenceCell< 0 > Vertex
int File_write_at_c(MPI_File fh, MPI_Offset offset, const void *buf, MPI_Count count, MPI_Datatype datatype, MPI_Status *status)
int File_write_at_all_c(MPI_File fh, MPI_Offset offset, const void *buf, MPI_Count count, MPI_Datatype datatype, MPI_Status *status)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
const MPI_Datatype mpi_type_id_for_type
void free_communicator(MPI_Comm mpi_communicator)
std::size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
unsigned int needed_digits(const unsigned int max_number)
constexpr T pow(const T base, const int iexp)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
constexpr types::global_dof_index invalid_size_type
constexpr unsigned int invalid_unsigned_int
::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 > &)
CFFlags(const double time=-std::numeric_limits< double >::infinity(), const bool keep_existing_file=false, const std::map< std::string, std::vector< std::pair< std::string, AttributeValue > > > &attributes={})
std::map< std::string, std::vector< std::pair< std::string, AttributeValue > > > attributes
std::variant< int, double, std::string > AttributeValue
void parse_parameters(const ParameterHandler &prm)
DXFlags(const bool write_neighbors=false, const bool int_binary=false, const bool coordinates_binary=false, const bool data_binary=false)
static void declare_parameters(ParameterHandler &prm)
static void declare_parameters(ParameterHandler &prm)
void parse_parameters(const ParameterHandler &prm)
bool filter_duplicate_vertices
DataOutFilterFlags(const bool filter_duplicate_vertices=false, const bool xdmf_hdf5_output=false)
static void declare_parameters(ParameterHandler &prm)
static RgbValues default_color_function(const double value, const double min_value, const double max_value)
void parse_parameters(const ParameterHandler &prm)
ColorFunction color_function
RgbValues(*)(const double value, const double min_value, const double max_value) ColorFunction
static RgbValues grey_scale_color_function(const double value, const double min_value, const double max_value)
EpsFlags(const unsigned int height_vector=0, const unsigned int color_vector=0, const SizeType size_type=width, const unsigned int size=300, const double line_width=0.5, const double azimut_angle=60, const double turn_angle=30, const double z_scaling=1.0, const bool draw_mesh=true, const bool draw_cells=true, const bool shade_cells=true, const ColorFunction color_function=&default_color_function)
unsigned int color_vector
static RgbValues reverse_grey_scale_color_function(const double value, const double min_value, const double max_value)
@ width
Scale to given width.
@ height
Scale to given height.
unsigned int height_vector
std::vector< std::string > space_dimension_labels
std::size_t memory_consumption() const
DataOutBase::CompressionLevel compression_level
Hdf5Flags(const CompressionLevel compression_level=CompressionLevel::best_speed)
static void declare_parameters(ParameterHandler &prm)
std::size_t memory_consumption() const
static const unsigned int no_neighbor
bool operator==(const Patch &patch) const
ReferenceCell< dim > reference_cell
void swap(Patch< dim, spacedim > &other_patch) noexcept
static const unsigned int space_dim
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 &prm)
PovrayFlags(const bool smooth=false, const bool bicubic_patch=false, const bool external_data=false)
void parse_parameters(const ParameterHandler &prm)
unsigned int height_vector
SvgFlags(const unsigned int height_vector=0, const int azimuth_angle=37, const int polar_angle=45, const unsigned int line_thickness=1, const bool margin=true, const bool draw_colorbar=true)
unsigned int line_thickness
std::size_t memory_consumption() const
TecplotFlags(const char *zone_name=nullptr, const double solution_time=-1.0)
void parse_parameters(const ParameterHandler &prm)
static void declare_parameters(ParameterHandler &prm)
UcdFlags(const bool write_preamble=false)
std::map< std::string, std::string > physical_units
bool write_higher_order_cells
DataOutBase::CompressionLevel compression_level
VtkFlags(const double time=std::numeric_limits< double >::lowest(), const unsigned int cycle=numbers::invalid_unsigned_int, const bool print_date_and_time=true, const CompressionLevel compression_level=CompressionLevel::best_speed, const bool write_higher_order_cells=false, const std::map< std::string, std::string > &physical_units={})
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices()
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)