40 #ifdef DEAL_II_WITH_ZLIB
44 #ifdef DEAL_II_WITH_HDF5
48 #include <boost/iostreams/copy.hpp>
49 #include <boost/iostreams/device/back_inserter.hpp>
50 #include <boost/iostreams/filtering_stream.hpp>
51 #ifdef DEAL_II_WITH_ZLIB
52 # include <boost/iostreams/filter/zlib.hpp>
66 <<
"Unexpected input: expected line\n <" << arg1
67 <<
">\nbut got\n <" << arg2 <<
">");
73 #ifdef DEAL_II_WITH_ZLIB
74 constexpr
bool deal_ii_with_zlib =
true;
76 constexpr
bool deal_ii_with_zlib =
false;
80 #ifdef DEAL_II_WITH_ZLIB
91 return Z_NO_COMPRESSION;
95 return Z_BEST_COMPRESSION;
97 return Z_DEFAULT_COMPRESSION;
100 return Z_NO_COMPRESSION;
104 # ifdef DEAL_II_WITH_MPI
115 return boost::iostreams::zlib::no_compression;
117 return boost::iostreams::zlib::best_speed;
119 return boost::iostreams::zlib::best_compression;
121 return boost::iostreams::zlib::default_compression;
124 return boost::iostreams::zlib::no_compression;
134 template <
typename T>
136 compress_array(
const std::vector<T> &data,
139 #ifdef DEAL_II_WITH_ZLIB
140 if (data.size() != 0)
142 const std::size_t uncompressed_size = (data.size() *
sizeof(
T));
154 auto compressed_data_length = compressBound(uncompressed_size);
159 std::vector<unsigned char> compressed_data(compressed_data_length);
161 int err = compress2(&compressed_data[0],
162 &compressed_data_length,
163 reinterpret_cast<const Bytef *
>(data.data()),
165 get_zlib_compression_level(compression_level));
170 compressed_data.resize(compressed_data_length);
173 const std::uint32_t compression_header[4] = {
175 static_cast<std::uint32_t
>(uncompressed_size),
176 static_cast<std::uint32_t
>(
178 static_cast<std::uint32_t
>(
179 compressed_data_length)};
181 const auto *
const header_start =
182 reinterpret_cast<const unsigned char *
>(&compression_header[0]);
185 {header_start, header_start + 4 *
sizeof(std::uint32_t)}) +
192 (void)compression_level;
194 ExcMessage(
"This function can only be called if cmake found "
195 "a working libz installation."));
210 template <
typename T>
212 vtu_stringize_array(
const std::vector<T> &data,
216 if (deal_ii_with_zlib &&
220 return compress_array(data, compression_level);
224 std::ostringstream stream;
225 stream.precision(precision);
226 for (
const T &el : data)
241 struct ParallelIntermediateHeader
244 std::uint64_t version;
245 std::uint64_t compression;
246 std::uint64_t dimension;
247 std::uint64_t space_dimension;
248 std::uint64_t n_ranks;
249 std::uint64_t n_patches;
357 template <
int dim,
int spacedim,
typename Number =
double>
358 std::unique_ptr<Table<2, Number>>
359 create_global_data_table(
const std::vector<Patch<dim, spacedim>> &patches)
363 return std::make_unique<Table<2, Number>>();
372 const unsigned int n_data_sets = patches[0].points_are_available ?
373 (patches[0].data.n_rows() - spacedim) :
374 patches[0].data.n_rows();
375 const unsigned int n_data_points =
376 std::accumulate(patches.begin(),
379 [](
const unsigned int count,
381 return count + patch.data.n_cols();
384 std::unique_ptr<Table<2, Number>> global_data_table =
385 std::make_unique<Table<2, Number>>(n_data_sets, n_data_points);
388 unsigned int next_value = 0;
389 for (
const auto &patch : patches)
391 const unsigned int n_subdivisions = patch.n_subdivisions;
392 (void)n_subdivisions;
394 Assert((patch.data.n_rows() == n_data_sets &&
395 !patch.points_are_available) ||
396 (patch.data.n_rows() == n_data_sets + spacedim &&
397 patch.points_are_available),
399 (n_data_sets + spacedim) :
401 patch.data.n_rows()));
402 Assert(patch.reference_cell != ReferenceCells::get_hypercube<dim>() ||
403 (n_data_sets == 0) ||
404 (patch.data.n_cols() ==
405 Utilities::fixed_power<dim>(n_subdivisions + 1)),
407 n_subdivisions + 1));
409 for (
unsigned int i = 0; i < patch.data.n_cols(); ++i, ++next_value)
410 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
411 (*global_data_table)[data_set][next_value] =
412 patch.data(data_set, i);
416 return global_data_table;
445 for (
unsigned int d = 0;
d < dim; ++
d)
449 unsigned int internal_ind;
460 internal_ind = it->second;
470 const unsigned int pt_index)
489 node_data[
node_dim * existing_point.second +
d] =
490 existing_point.first(
d);
498 std::vector<unsigned int> &cell_data)
const
504 cell_data[filtered_cell.first] =
505 filtered_cell.second + local_node_offset;
574 const unsigned int start,
575 const std::array<unsigned int, dim> &offsets)
579 const unsigned int base_entry =
592 const unsigned int d1 = offsets[0];
601 const unsigned int d1 = offsets[0];
602 const unsigned int d2 = offsets[1];
613 const unsigned int d1 = offsets[0];
614 const unsigned int d2 = offsets[1];
615 const unsigned int d3 = offsets[2];
637 const unsigned int start,
638 const unsigned int n_points,
643 const unsigned int base_entry = index * n_points;
645 static const std::array<unsigned int, 5> table = {{0, 1, 3, 2, 4}};
647 for (
unsigned int i = 0; i < n_points; ++i)
658 const unsigned int dimension,
659 const unsigned int set_num,
662 unsigned int new_dim;
681 for (
unsigned int d = 0;
d < new_dim; ++
d)
684 data_sets.back()[r * new_dim +
d] = data_vectors(set_num +
d, i);
700 const char *gmv_cell_type[4] = {
"",
"line 2",
"quad 4",
"hex 8"};
702 const char *ucd_cell_type[4] = {
"pt",
"line",
"quad",
"hex"};
704 const char *tecplot_cell_type[4] = {
"",
"lineseg",
"quadrilateral",
"brick"};
719 template <
int dim,
int spacedim>
720 std::array<unsigned int, 3>
722 const bool write_higher_order_cells)
724 std::array<unsigned int, 3> vtk_cell_id = {
730 if (write_higher_order_cells)
733 vtk_cell_id[2] = patch.
data.n_cols();
739 patch.
data.n_cols() == 6)
743 vtk_cell_id[2] = patch.
data.n_cols();
746 patch.
data.n_cols() == 10)
750 vtk_cell_id[2] = patch.
data.n_cols();
775 template <
int dim,
int spacedim>
777 get_equispaced_location(
779 const std::initializer_list<unsigned int> &lattice_location,
780 const unsigned int n_subdivisions)
787 const unsigned int xstep = (dim > 0 ? *(lattice_location.begin() + 0) : 0);
788 const unsigned int ystep = (dim > 1 ? *(lattice_location.begin() + 1) : 0);
789 const unsigned int zstep = (dim > 2 ? *(lattice_location.begin() + 2) : 0);
797 unsigned int point_no = 0;
802 point_no += (n_subdivisions + 1) * (n_subdivisions + 1) * zstep;
806 point_no += (n_subdivisions + 1) * ystep;
820 for (
unsigned int d = 0;
d < spacedim; ++
d)
821 node[
d] = patch.
data(patch.
data.size(0) - spacedim +
d, point_no);
833 const double stepsize = 1. / n_subdivisions;
834 const double xfrac = xstep * stepsize;
840 const double yfrac = ystep * stepsize;
842 node += ((patch.
vertices[3] * xfrac) +
843 (patch.
vertices[2] * (1 - xfrac))) *
847 const double zfrac = zstep * stepsize;
849 node += (((patch.
vertices[5] * xfrac) +
850 (patch.
vertices[4] * (1 - xfrac))) *
853 (patch.
vertices[6] * (1 - xfrac))) *
866 template <
int dim,
int spacedim>
869 const unsigned int node_index)
874 unsigned int point_no_actual = node_index;
879 static const std::array<unsigned int, 5> table = {{0, 1, 3, 2, 4}};
880 point_no_actual = table[node_index];
888 for (
unsigned int d = 0;
d < spacedim; ++
d)
890 patch.
data(patch.
data.size(0) - spacedim +
d, point_no_actual);
903 return patch.
vertices[point_no_actual];
914 template <
int dim,
int spacedim>
915 std::tuple<unsigned int, unsigned int>
916 count_nodes_and_cells(
919 unsigned int n_nodes = 0;
921 for (
const auto &patch : patches)
925 "The reference cell for this patch is set to 'Invalid', "
926 "but that is clearly not a valid choice. Did you forget "
927 "to set the reference cell for the patch?"));
942 return std::make_tuple(n_nodes,
n_cells);
952 template <
int dim,
int spacedim>
953 std::tuple<unsigned int, unsigned int, unsigned int>
954 count_nodes_and_cells_and_points(
956 const bool write_higher_order_cells)
958 unsigned int n_nodes = 0;
960 unsigned int n_points_and_n_cells = 0;
962 for (
const auto &patch : patches)
968 if (write_higher_order_cells)
974 n_points_and_n_cells +=
983 const unsigned int n_subcells =
986 n_points_and_n_cells +=
992 n_nodes += patch.
data.n_cols();
994 n_points_and_n_cells += patch.
data.n_cols() + 1;
998 return std::make_tuple(n_nodes,
n_cells, n_points_and_n_cells);
1006 template <
typename FlagsType>
1013 StreamBase(std::ostream &stream,
const FlagsType &flags)
1025 write_point(
const unsigned int,
const Point<dim> &)
1028 ExcMessage(
"The derived class you are using needs to "
1029 "reimplement this function if you want to call "
1049 write_cell(
const unsigned int ,
1050 const unsigned int ,
1051 std::array<unsigned int, dim> & )
1054 ExcMessage(
"The derived class you are using needs to "
1055 "reimplement this function if you want to call "
1066 write_cell_single(
const unsigned int index,
1067 const unsigned int start,
1068 const unsigned int n_points,
1077 ExcMessage(
"The derived class you are using needs to "
1078 "reimplement this function if you want to call "
1096 template <
typename T>
1110 unsigned int selected_component;
1117 std::ostream &stream;
1122 const FlagsType flags;
1128 class DXStream :
public StreamBase<DataOutBase::DXFlags>
1135 write_point(
const unsigned int index,
const Point<dim> &);
1147 write_cell(
const unsigned int index,
1148 const unsigned int start,
1149 const std::array<unsigned int, dim> &offsets);
1157 template <
typename data>
1159 write_dataset(
const unsigned int index,
const std::vector<data> &
values);
1165 class GmvStream :
public StreamBase<DataOutBase::GmvFlags>
1172 write_point(
const unsigned int index,
const Point<dim> &);
1184 write_cell(
const unsigned int index,
1185 const unsigned int start,
1186 const std::array<unsigned int, dim> &offsets);
1192 class TecplotStream :
public StreamBase<DataOutBase::TecplotFlags>
1199 write_point(
const unsigned int index,
const Point<dim> &);
1211 write_cell(
const unsigned int index,
1212 const unsigned int start,
1213 const std::array<unsigned int, dim> &offsets);
1219 class UcdStream :
public StreamBase<DataOutBase::UcdFlags>
1226 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);
1250 template <
typename data>
1252 write_dataset(
const unsigned int index,
const std::vector<data> &
values);
1258 class VtkStream :
public StreamBase<DataOutBase::VtkFlags>
1265 write_point(
const unsigned int index,
const Point<dim> &);
1277 write_cell(
const unsigned int index,
1278 const unsigned int start,
1279 const std::array<unsigned int, dim> &offsets);
1285 write_cell_single(
const unsigned int index,
1286 const unsigned int start,
1287 const unsigned int n_points,
1299 write_high_order_cell(
const unsigned int start,
1300 const std::vector<unsigned> &connectivity);
1313 DXStream::write_point(
const unsigned int,
const Point<dim> &p)
1315 if (flags.coordinates_binary)
1318 for (
unsigned int d = 0;
d < dim; ++
d)
1320 stream.write(
reinterpret_cast<const char *
>(data), dim *
sizeof(*data));
1324 for (
unsigned int d = 0;
d < dim; ++
d)
1325 stream << p(
d) <<
'\t';
1339 std::array<unsigned int, GeometryInfo<0>::vertices_per_cell>
1340 set_node_numbers(
const unsigned int ,
1341 const std::array<unsigned int, 0> & )
1349 std::array<unsigned int, GeometryInfo<1>::vertices_per_cell>
1350 set_node_numbers(
const unsigned int start,
1351 const std::array<unsigned int, 1> &offsets)
1353 std::array<unsigned int, GeometryInfo<1>::vertices_per_cell> nodes;
1355 nodes[1] = start + offsets[0];
1361 std::array<unsigned int, GeometryInfo<2>::vertices_per_cell>
1362 set_node_numbers(
const unsigned int start,
1363 const std::array<unsigned int, 2> &offsets)
1366 const unsigned int d1 = offsets[0];
1367 const unsigned int d2 = offsets[1];
1369 std::array<unsigned int, GeometryInfo<2>::vertices_per_cell> nodes;
1371 nodes[1] = start + d1;
1372 nodes[2] = start + d2;
1373 nodes[3] = start + d2 + d1;
1379 std::array<unsigned int, GeometryInfo<3>::vertices_per_cell>
1380 set_node_numbers(
const unsigned int start,
1381 const std::array<unsigned int, 3> &offsets)
1383 const unsigned int d1 = offsets[0];
1384 const unsigned int d2 = offsets[1];
1385 const unsigned int d3 = offsets[2];
1387 std::array<unsigned int, GeometryInfo<3>::vertices_per_cell> nodes;
1389 nodes[1] = start + d1;
1390 nodes[2] = start + d2;
1391 nodes[3] = start + d2 + d1;
1392 nodes[4] = start + d3;
1393 nodes[5] = start + d3 + d1;
1394 nodes[6] = start + d3 + d2;
1395 nodes[7] = start + d3 + d2 + d1;
1404 DXStream::write_cell(
const unsigned int,
1405 const unsigned int start,
1406 const std::array<unsigned int, dim> &offsets)
1409 DataOutBaseImplementation::set_node_numbers(start, offsets);
1411 if (flags.int_binary)
1413 std::array<unsigned int, GeometryInfo<dim>::vertices_per_cell> temp;
1414 for (
unsigned int i = 0; i < nodes.size(); ++i)
1416 stream.write(
reinterpret_cast<const char *
>(temp.data()),
1417 temp.size() *
sizeof(temp[0]));
1421 for (
unsigned int i = 0; i < nodes.size() - 1; ++i)
1423 stream << nodes[GeometryInfo<dim>::dx_to_deal[nodes.size() - 1]]
1430 template <
typename data>
1432 DXStream::write_dataset(
const unsigned int,
const std::vector<data> &
values)
1434 if (flags.data_binary)
1436 stream.write(
reinterpret_cast<const char *
>(
values.data()),
1437 values.size() *
sizeof(data));
1441 for (
unsigned int i = 0; i <
values.size(); ++i)
1442 stream <<
'\t' <<
values[i];
1458 GmvStream::write_point(
const unsigned int,
const Point<dim> &p)
1462 stream << p(selected_component) <<
' ';
1469 GmvStream::write_cell(
const unsigned int,
1470 const unsigned int s,
1471 const std::array<unsigned int, dim> &offsets)
1474 const unsigned int start = s + 1;
1475 stream << gmv_cell_type[dim] <<
'\n';
1487 const unsigned int d1 = offsets[0];
1489 stream <<
'\t' << start + d1;
1495 const unsigned int d1 = offsets[0];
1496 const unsigned int d2 = offsets[1];
1498 stream <<
'\t' << start + d1;
1499 stream <<
'\t' << start + d2 + d1 <<
'\t' << start + d2;
1505 const unsigned int d1 = offsets[0];
1506 const unsigned int d2 = offsets[1];
1507 const unsigned int d3 = offsets[2];
1509 stream <<
'\t' << start + d1;
1510 stream <<
'\t' << start + d2 + d1 <<
'\t' << start + d2;
1511 stream <<
'\t' << start + d3 <<
'\t' << start + d3 + d1 <<
'\t'
1512 << start + d3 + d2 + d1 <<
'\t' << start + d3 + d2;
1524 TecplotStream::TecplotStream(std::ostream &out,
1532 TecplotStream::write_point(
const unsigned int,
const Point<dim> &p)
1536 stream << p(selected_component) <<
'\n';
1543 TecplotStream::write_cell(
const unsigned int,
1544 const unsigned int s,
1545 const std::array<unsigned int, dim> &offsets)
1547 const unsigned int start = s + 1;
1559 const unsigned int d1 = offsets[0];
1561 stream <<
'\t' << start + d1;
1567 const unsigned int d1 = offsets[0];
1568 const unsigned int d2 = offsets[1];
1570 stream <<
'\t' << start + d1;
1571 stream <<
'\t' << start + d2 + d1 <<
'\t' << start + d2;
1577 const unsigned int d1 = offsets[0];
1578 const unsigned int d2 = offsets[1];
1579 const unsigned int d3 = offsets[2];
1581 stream <<
'\t' << start + d1;
1582 stream <<
'\t' << start + d2 + d1 <<
'\t' << start + d2;
1583 stream <<
'\t' << start + d3 <<
'\t' << start + d3 + d1 <<
'\t'
1584 << start + d3 + d2 + d1 <<
'\t' << start + d3 + d2;
1603 UcdStream::write_point(
const unsigned int index,
const Point<dim> &p)
1605 stream <<
index + 1 <<
" ";
1607 for (
unsigned int i = 0; i < dim; ++i)
1608 stream << p(i) <<
' ';
1610 for (
unsigned int i = dim; i < 3; ++i)
1619 UcdStream::write_cell(
const unsigned int index,
1620 const unsigned int start,
1621 const std::array<unsigned int, dim> &offsets)
1624 DataOutBaseImplementation::set_node_numbers(start, offsets);
1627 stream <<
index + 1 <<
"\t0 " << ucd_cell_type[dim];
1628 for (
unsigned int i = 0; i < nodes.size(); ++i)
1635 template <
typename data>
1637 UcdStream::write_dataset(
const unsigned int index,
1638 const std::vector<data> &
values)
1640 stream <<
index + 1;
1641 for (
unsigned int i = 0; i <
values.size(); ++i)
1642 stream <<
'\t' <<
values[i];
1657 VtkStream::write_point(
const unsigned int,
const Point<dim> &p)
1662 for (
unsigned int i = dim; i < 3; ++i)
1671 VtkStream::write_cell(
const unsigned int,
1672 const unsigned int start,
1673 const std::array<unsigned int, dim> &offsets)
1675 stream << GeometryInfo<dim>::vertices_per_cell <<
'\t';
1687 const unsigned int d1 = offsets[0];
1689 stream <<
'\t' << start + d1;
1695 const unsigned int d1 = offsets[0];
1696 const unsigned int d2 = offsets[1];
1698 stream <<
'\t' << start + d1;
1699 stream <<
'\t' << start + d2 + d1 <<
'\t' << start + d2;
1705 const unsigned int d1 = offsets[0];
1706 const unsigned int d2 = offsets[1];
1707 const unsigned int d3 = offsets[2];
1709 stream <<
'\t' << start + d1;
1710 stream <<
'\t' << start + d2 + d1 <<
'\t' << start + d2;
1711 stream <<
'\t' << start + d3 <<
'\t' << start + d3 + d1 <<
'\t'
1712 << start + d3 + d2 + d1 <<
'\t' << start + d3 + d2;
1725 VtkStream::write_cell_single(
const unsigned int index,
1726 const unsigned int start,
1727 const unsigned int n_points,
1732 static const std::array<unsigned int, 5> table = {{0, 1, 3, 2, 4}};
1734 stream <<
'\t' << n_points;
1735 for (
unsigned int i = 0; i < n_points; ++i)
1744 VtkStream::write_high_order_cell(
const unsigned int start,
1745 const std::vector<unsigned> &connectivity)
1747 stream << connectivity.size();
1748 for (
const auto &c : connectivity)
1749 stream <<
'\t' << start + c;
1761 template <
int dim,
int spacedim>
1765 template <
int dim,
int spacedim>
1769 template <
int dim,
int spacedim>
1771 : patch_index(no_neighbor)
1773 , points_are_available(false)
1787 template <
int dim,
int spacedim>
1813 if (data.n_rows() != patch.
data.n_rows())
1816 if (data.n_cols() != patch.
data.n_cols())
1819 for (
unsigned int i = 0; i < data.n_rows(); ++i)
1820 for (
unsigned int j = 0; j < data.n_cols(); ++j)
1821 if (data[i][j] != patch.
data[i][j])
1829 template <
int dim,
int spacedim>
1835 sizeof(neighbors) /
sizeof(neighbors[0]) *
1846 template <
int dim,
int spacedim>
1854 data.swap(other_patch.
data);
1861 template <
int spacedim>
1865 template <
int spacedim>
1869 template <
int spacedim>
1873 template <
int spacedim>
1876 template <
int spacedim>
1880 template <
int spacedim>
1882 : patch_index(no_neighbor)
1883 , points_are_available(false)
1890 template <
int spacedim>
1894 const unsigned int dim = 0;
1908 if (
data.n_rows() != patch.
data.n_rows())
1911 if (
data.n_cols() != patch.
data.n_cols())
1914 for (
unsigned int i = 0; i <
data.n_rows(); ++i)
1915 for (
unsigned int j = 0; j <
data.n_cols(); ++j)
1916 if (
data[i][j] != patch.
data[i][j])
1924 template <
int spacedim>
1936 template <
int spacedim>
1949 : write_preamble(write_preamble)
1964 : space_dimension_labels(labels)
1978 const bool bicubic_patch,
1979 const bool external_data)
1981 , bicubic_patch(bicubic_patch)
1982 , external_data(external_data)
1987 const bool xdmf_hdf5_output)
1988 : filter_duplicate_vertices(filter_duplicate_vertices)
1989 , xdmf_hdf5_output(xdmf_hdf5_output)
1997 "Filter duplicate vertices",
2000 "Whether to remove duplicate vertex values. deal.II duplicates "
2001 "vertices once for each adjacent cell so that it can output "
2002 "discontinuous quantities for which there may be more than one "
2003 "value for each vertex position. Setting this flag to "
2004 "'true' will merge all of these values by selecting a "
2005 "random one and outputting this as 'the' value for the vertex. "
2006 "As long as the data to be output corresponds to continuous "
2007 "fields, merging vertices has no effect. On the other hand, "
2008 "if the data to be output corresponds to discontinuous fields "
2009 "(either because you are using a discontinuous finite element, "
2010 "or because you are using a DataPostprocessor that yields "
2011 "discontinuous data, or because the data to be output has been "
2012 "produced by entirely different means), then the data in the "
2013 "output file no longer faithfully represents the underlying data "
2014 "because the discontinuous field has been replaced by a "
2015 "continuous one. Note also that the filtering can not occur "
2016 "on processor boundaries. Thus, a filtered discontinuous field "
2017 "looks like a continuous field inside of a subdomain, "
2018 "but like a discontinuous field at the subdomain boundary."
2020 "In any case, filtering results in drastically smaller output "
2021 "files (smaller by about a factor of 2^dim).");
2026 "Whether the data will be used in an XDMF/HDF5 combination.");
2041 const bool int_binary,
2042 const bool coordinates_binary,
2043 const bool data_binary)
2044 : write_neighbors(write_neighbors)
2045 , int_binary(int_binary)
2046 , coordinates_binary(coordinates_binary)
2047 , data_binary(data_binary)
2048 , data_double(false)
2058 "A boolean field indicating whether neighborship "
2059 "information between cells is to be written to the "
2060 "OpenDX output file");
2064 "Output format of integer numbers, which is "
2065 "either a text representation (ascii) or binary integer "
2066 "values of 32 or 64 bits length");
2070 "Output format of vertex coordinates, which is "
2071 "either a text representation (ascii) or binary "
2072 "floating point values of 32 or 64 bits length");
2076 "Output format of data values, which is "
2077 "either a text representation (ascii) or binary "
2078 "floating point values of 32 or 64 bits length");
2098 "A flag indicating whether a comment should be "
2099 "written to the beginning of the output file "
2100 "indicating date and time of creation as well "
2101 "as the creating program");
2115 const int azimuth_angle,
2116 const int polar_angle,
2117 const unsigned int line_thickness,
2119 const bool draw_colorbar)
2122 , height_vector(height_vector)
2123 , azimuth_angle(azimuth_angle)
2124 , polar_angle(polar_angle)
2125 , line_thickness(line_thickness)
2127 , draw_colorbar(draw_colorbar)
2138 "A flag indicating whether POVRAY should use smoothed "
2139 "triangles instead of the usual ones");
2143 "Whether POVRAY should use bicubic patches");
2147 "Whether camera and lighting information should "
2148 "be put into an external file \"data.inc\" or into "
2149 "the POVRAY input file");
2165 const unsigned int color_vector,
2167 const unsigned int size,
2168 const double line_width,
2169 const double azimut_angle,
2170 const double turn_angle,
2171 const double z_scaling,
2172 const bool draw_mesh,
2173 const bool draw_cells,
2174 const bool shade_cells,
2176 : height_vector(height_vector)
2177 , color_vector(color_vector)
2180 , line_width(line_width)
2181 , azimut_angle(azimut_angle)
2182 , turn_angle(turn_angle)
2183 , z_scaling(z_scaling)
2184 , draw_mesh(draw_mesh)
2185 , draw_cells(draw_cells)
2186 , shade_cells(shade_cells)
2187 , color_function(color_function)
2226 double sum = xmax + xmin;
2227 double sum13 = xmin + 3 * xmax;
2228 double sum22 = 2 * xmin + 2 * xmax;
2229 double sum31 = 3 * xmin + xmax;
2230 double dif = xmax - xmin;
2231 double rezdif = 1.0 / dif;
2235 if (x < (sum31) / 4)
2237 else if (x < (sum22) / 4)
2239 else if (x < (sum13) / 4)
2250 rgb_values.
green = 0;
2251 rgb_values.
blue = (x - xmin) * 4. * rezdif;
2255 rgb_values.
green = (4 * x - 3 * xmin - xmax) * rezdif;
2256 rgb_values.
blue = (sum22 - 4. * x) * rezdif;
2259 rgb_values.
red = (4 * x - 2 *
sum) * rezdif;
2260 rgb_values.
green = (xmin + 3 * xmax - 4 * x) * rezdif;
2261 rgb_values.
blue = 0;
2265 rgb_values.
green = (4 * x - xmin - 3 * xmax) * rezdif;
2266 rgb_values.
blue = (4. * x - sum13) * rezdif;
2273 rgb_values.
red = rgb_values.
green = rgb_values.
blue = 1;
2287 (x - xmin) / (xmax - xmin);
2300 1 - (x - xmin) / (xmax - xmin);
2312 "Number of the input vector that is to be used to "
2313 "generate height information");
2317 "Number of the input vector that is to be used to "
2318 "generate color information");
2322 "Whether width or height should be scaled to match "
2327 "The size (width or height) to which the eps output "
2328 "file is to be scaled");
2332 "The width in which the postscript renderer is to "
2337 "Angle of the viewing position against the vertical "
2342 "Angle of the viewing direction against the y-axis");
2346 "Scaling for the z-direction relative to the scaling "
2347 "used in x- and y-directions");
2351 "Whether the mesh lines, or only the surface should be "
2356 "Whether only the mesh lines, or also the interior of "
2357 "cells should be plotted. If this flag is false, then "
2358 "one can see through the mesh");
2362 "Whether the interior of cells shall be shaded");
2366 "default|grey scale|reverse grey scale"),
2367 "Name of a color function used to colorize mesh lines "
2368 "and/or cell interiors");
2378 if (prm.
get(
"Scale to width or height") ==
"width")
2390 if (prm.
get(
"Color function") ==
"default")
2392 else if (prm.
get(
"Color function") ==
"grey scale")
2394 else if (prm.
get(
"Color function") ==
"reverse grey scale")
2404 : compression_level(compression_level)
2409 : zone_name(zone_name)
2410 , solution_time(solution_time)
2424 const unsigned int cycle,
2425 const bool print_date_and_time,
2427 const bool write_higher_order_cells,
2428 const std::map<std::string, std::string> &physical_units)
2431 , print_date_and_time(print_date_and_time)
2432 , compression_level(compression_level)
2433 , write_higher_order_cells(write_higher_order_cells)
2434 , physical_units(physical_units)
2442 if (format_name ==
"none")
2445 if (format_name ==
"dx")
2448 if (format_name ==
"ucd")
2451 if (format_name ==
"gnuplot")
2454 if (format_name ==
"povray")
2457 if (format_name ==
"eps")
2460 if (format_name ==
"gmv")
2463 if (format_name ==
"tecplot")
2466 if (format_name ==
"vtk")
2469 if (format_name ==
"vtu")
2472 if (format_name ==
"deal.II intermediate")
2475 if (format_name ==
"hdf5")
2479 ExcMessage(
"The given file format name is not recognized: <" +
2480 format_name +
">"));
2491 return "none|dx|ucd|gnuplot|povray|eps|gmv|tecplot|vtk|vtu|hdf5|svg|deal.II intermediate";
2499 switch (output_format)
2541 template <
int dim,
int spacedim>
2542 std::vector<Point<spacedim>>
2546 static const std::array<unsigned int, 5> table = {{0, 1, 3, 2, 4}};
2548 std::vector<Point<spacedim>> node_positions;
2549 for (
const auto &patch : patches)
2552 if (patch.
reference_cell != ReferenceCells::get_hypercube<dim>())
2554 for (
unsigned int point_no = 0; point_no < patch.
data.n_cols();
2556 node_positions.emplace_back(get_node_location(
2565 const unsigned int n = n_subdivisions + 1;
2570 node_positions.emplace_back(
2571 get_equispaced_location(patch, {}, n_subdivisions));
2574 for (
unsigned int i1 = 0; i1 < n; ++i1)
2575 node_positions.emplace_back(
2576 get_equispaced_location(patch, {i1}, n_subdivisions));
2579 for (
unsigned int i2 = 0; i2 < n; ++i2)
2580 for (
unsigned int i1 = 0; i1 < n; ++i1)
2581 node_positions.emplace_back(get_equispaced_location(
2582 patch, {i1, i2}, n_subdivisions));
2585 for (
unsigned int i3 = 0; i3 < n; ++i3)
2586 for (
unsigned int i2 = 0; i2 < n; ++i2)
2587 for (
unsigned int i1 = 0; i1 < n; ++i1)
2588 node_positions.emplace_back(get_equispaced_location(
2589 patch, {i1, i2, i3}, n_subdivisions));
2598 return node_positions;
2602 template <
int dim,
int spacedim,
typename StreamType>
2608 const std::vector<Point<spacedim>> node_positions =
2612 for (
const auto &node : node_positions)
2613 out.write_point(count++, node);
2619 template <
int dim,
int spacedim,
typename StreamType>
2624 unsigned int count = 0;
2625 unsigned int first_vertex_of_patch = 0;
2626 for (
const auto &patch : patches)
2629 if (patch.
reference_cell != ReferenceCells::get_hypercube<dim>())
2631 out.write_cell_single(count++,
2632 first_vertex_of_patch,
2633 patch.
data.n_cols(),
2635 first_vertex_of_patch += patch.
data.n_cols();
2640 const unsigned int n = n_subdivisions + 1;
2646 const unsigned int offset = first_vertex_of_patch;
2647 out.template write_cell<0>(count++, offset, {});
2653 constexpr
unsigned int d1 = 1;
2655 for (
unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
2657 const unsigned int offset =
2658 first_vertex_of_patch + i1 * d1;
2659 out.template write_cell<1>(count++, offset, {{d1}});
2667 constexpr
unsigned int d1 = 1;
2668 const unsigned int d2 = n;
2670 for (
unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
2671 for (
unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
2673 const unsigned int offset =
2674 first_vertex_of_patch + i2 * d2 + i1 * d1;
2675 out.template write_cell<2>(count++,
2685 constexpr
unsigned int d1 = 1;
2686 const unsigned int d2 = n;
2687 const unsigned int d3 = n * n;
2689 for (
unsigned int i3 = 0; i3 < n_subdivisions; ++i3)
2690 for (
unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
2691 for (
unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
2693 const unsigned int offset = first_vertex_of_patch +
2696 out.template write_cell<3>(count++,
2708 first_vertex_of_patch +=
2709 Utilities::fixed_power<dim>(n_subdivisions + 1);
2718 template <
int dim,
int spacedim,
typename StreamType>
2722 const bool legacy_format)
2725 unsigned int first_vertex_of_patch = 0;
2727 std::vector<unsigned> connectivity;
2729 for (
const auto &patch : patches)
2731 if (patch.
reference_cell != ReferenceCells::get_hypercube<dim>())
2733 connectivity.resize(patch.
data.n_cols());
2735 for (
unsigned int i = 0; i < patch.
data.n_cols(); ++i)
2736 connectivity[i] = i;
2738 out.template write_high_order_cell<dim>(first_vertex_of_patch,
2741 first_vertex_of_patch += patch.
data.n_cols();
2746 const unsigned int n = n_subdivisions + 1;
2748 connectivity.resize(Utilities::fixed_power<dim>(n));
2755 ExcMessage(
"Point-like cells should not be possible "
2756 "when writing higher-order cells."));
2761 for (
unsigned int i1 = 0; i1 < n_subdivisions + 1; ++i1)
2763 const unsigned int local_index = i1;
2764 const unsigned int connectivity_index =
2766 .template vtk_lexicographic_to_node_index<1>(
2767 {{i1}}, {{n_subdivisions}}, legacy_format);
2768 connectivity[connectivity_index] = local_index;
2775 for (
unsigned int i2 = 0; i2 < n_subdivisions + 1; ++i2)
2776 for (
unsigned int i1 = 0; i1 < n_subdivisions + 1; ++i1)
2778 const unsigned int local_index = i2 * n + i1;
2779 const unsigned int connectivity_index =
2781 .template vtk_lexicographic_to_node_index<2>(
2783 {{n_subdivisions, n_subdivisions}},
2785 connectivity[connectivity_index] = local_index;
2792 for (
unsigned int i3 = 0; i3 < n_subdivisions + 1; ++i3)
2793 for (
unsigned int i2 = 0; i2 < n_subdivisions + 1; ++i2)
2794 for (
unsigned int i1 = 0; i1 < n_subdivisions + 1; ++i1)
2796 const unsigned int local_index =
2797 i3 * n * n + i2 * n + i1;
2798 const unsigned int connectivity_index =
2800 .template vtk_lexicographic_to_node_index<3>(
2806 connectivity[connectivity_index] = local_index;
2817 out.template write_high_order_cell<dim>(first_vertex_of_patch,
2821 first_vertex_of_patch += Utilities::fixed_power<dim>(n);
2829 template <
int dim,
int spacedim,
typename StreamType>
2832 unsigned int n_data_sets,
2833 const bool double_precision,
2837 unsigned int count = 0;
2839 for (
const auto &patch : patches)
2842 const unsigned int n = n_subdivisions + 1;
2846 (patch.
data.n_rows() == n_data_sets + spacedim &&
2849 (n_data_sets + spacedim) :
2851 patch.
data.n_rows()));
2852 Assert(patch.
data.n_cols() == Utilities::fixed_power<dim>(n),
2855 std::vector<float> floats(n_data_sets);
2856 std::vector<double> doubles(n_data_sets);
2859 for (
unsigned int i = 0; i < Utilities::fixed_power<dim>(n);
2861 if (double_precision)
2863 for (
unsigned int data_set = 0; data_set < n_data_sets;
2865 doubles[data_set] = patch.
data(data_set, i);
2866 out.write_dataset(count, doubles);
2870 for (
unsigned int data_set = 0; data_set < n_data_sets;
2872 floats[data_set] = patch.
data(data_set, i);
2873 out.write_dataset(count, floats);
2898 camera_vertical[0] = camera_horizontal[1] * camera_direction[2] -
2899 camera_horizontal[2] * camera_direction[1];
2900 camera_vertical[1] = camera_horizontal[2] * camera_direction[0] -
2901 camera_horizontal[0] * camera_direction[2];
2902 camera_vertical[2] = camera_horizontal[0] * camera_direction[1] -
2903 camera_horizontal[1] * camera_direction[0];
2907 phi /= (
point[0] - camera_position[0]) * camera_direction[0] +
2908 (
point[1] - camera_position[1]) * camera_direction[1] +
2909 (
point[2] - camera_position[2]) * camera_direction[2];
2913 camera_position[0] + phi * (
point[0] - camera_position[0]);
2915 camera_position[1] + phi * (
point[1] - camera_position[1]);
2917 camera_position[2] + phi * (
point[2] - camera_position[2]);
2920 projection_decomposition[0] = (projection[0] - camera_position[0] -
2921 camera_focus * camera_direction[0]) *
2922 camera_horizontal[0];
2923 projection_decomposition[0] += (projection[1] - camera_position[1] -
2924 camera_focus * camera_direction[1]) *
2925 camera_horizontal[1];
2926 projection_decomposition[0] += (projection[2] - camera_position[2] -
2927 camera_focus * camera_direction[2]) *
2928 camera_horizontal[2];
2930 projection_decomposition[1] = (projection[0] - camera_position[0] -
2931 camera_focus * camera_direction[0]) *
2933 projection_decomposition[1] += (projection[1] - camera_position[1] -
2934 camera_focus * camera_direction[1]) *
2936 projection_decomposition[1] += (projection[2] - camera_position[2] -
2937 camera_focus * camera_direction[2]) *
2940 return projection_decomposition;
2949 svg_get_gradient_parameters(
Point<3> points[])
2955 for (
int i = 0; i < 2; ++i)
2957 for (
int j = 0; j < 2 - i; ++j)
2959 if (points[j][2] > points[j + 1][2])
2962 points[j] = points[j + 1];
2963 points[j + 1] = temp;
2970 v_inter = points[1];
2977 A[0][0] = v_max[0] - v_min[0];
2978 A[0][1] = v_inter[0] - v_min[0];
2979 A[1][0] = v_max[1] - v_min[1];
2980 A[1][1] = v_inter[1] - v_min[1];
2986 bool col_change =
false;
2995 double temp =
A[1][0];
3000 for (
unsigned int k = 0; k < 1; ++k)
3002 for (
unsigned int i = k + 1; i < 2; ++i)
3004 x =
A[i][k] /
A[k][k];
3006 for (
unsigned int j = k + 1; j < 2; ++j)
3007 A[i][j] =
A[i][j] -
A[k][j] * x;
3009 b[i] =
b[i] -
b[k] * x;
3013 b[1] =
b[1] /
A[1][1];
3015 for (
int i = 0; i >= 0; i--)
3019 for (
unsigned int j = i + 1; j < 2; ++j)
3022 b[i] =
sum /
A[i][i];
3032 double c =
b[0] * (v_max[2] - v_min[2]) +
b[1] * (v_inter[2] - v_min[2]) +
3036 A[0][0] = v_max[0] - v_min[0];
3037 A[0][1] = v_inter[0] - v_min[0];
3038 A[1][0] = v_max[1] - v_min[1];
3039 A[1][1] = v_inter[1] - v_min[1];
3041 b[0] = 1.0 - v_min[0];
3053 double temp =
A[1][0];
3058 for (
unsigned int k = 0; k < 1; ++k)
3060 for (
unsigned int i = k + 1; i < 2; ++i)
3062 x =
A[i][k] /
A[k][k];
3064 for (
unsigned int j = k + 1; j < 2; ++j)
3065 A[i][j] =
A[i][j] -
A[k][j] * x;
3067 b[i] =
b[i] -
b[k] * x;
3071 b[1] =
b[1] /
A[1][1];
3073 for (
int i = 0; i >= 0; i--)
3077 for (
unsigned int j = i + 1; j < 2; ++j)
3080 b[i] =
sum /
A[i][i];
3090 gradient[0] =
b[0] * (v_max[2] - v_min[2]) +
3091 b[1] * (v_inter[2] - v_min[2]) - c + v_min[2];
3094 A[0][0] = v_max[0] - v_min[0];
3095 A[0][1] = v_inter[0] - v_min[0];
3096 A[1][0] = v_max[1] - v_min[1];
3097 A[1][1] = v_inter[1] - v_min[1];
3100 b[1] = 1.0 - v_min[1];
3111 double temp =
A[1][0];
3116 for (
unsigned int k = 0; k < 1; ++k)
3118 for (
unsigned int i = k + 1; i < 2; ++i)
3120 x =
A[i][k] /
A[k][k];
3122 for (
unsigned int j = k + 1; j < 2; ++j)
3123 A[i][j] =
A[i][j] -
A[k][j] * x;
3125 b[i] =
b[i] -
b[k] * x;
3129 b[1] =
b[1] /
A[1][1];
3131 for (
int i = 0; i >= 0; i--)
3135 for (
unsigned int j = i + 1; j < 2; ++j)
3138 b[i] =
sum /
A[i][i];
3148 gradient[1] =
b[0] * (v_max[2] - v_min[2]) +
3149 b[1] * (v_inter[2] - v_min[2]) - c + v_min[2];
3155 gradient[1] * (v_min[1] - v_max[1]);
3159 gradient_parameters[0] = v_min[0];
3160 gradient_parameters[1] = v_min[1];
3165 gradient_parameters[4] = v_min[2];
3166 gradient_parameters[5] = v_max[2];
3168 return gradient_parameters;
3174 template <
int dim,
int spacedim>
3178 const std::vector<std::string> &data_names,
3180 std::tuple<
unsigned int,
3193 #ifndef DEAL_II_WITH_MPI
3202 if (patches.empty())
3206 const unsigned int n_data_sets = data_names.size();
3208 UcdStream ucd_out(out, flags);
3211 unsigned int n_nodes;
3213 std::tie(n_nodes,
n_cells) = count_nodes_and_cells(patches);
3219 <<
"# This file was generated by the deal.II library." <<
'\n'
3223 <<
"# For a description of the UCD format see the AVS Developer's guide."
3229 out << n_nodes <<
' ' <<
n_cells <<
' ' << n_data_sets <<
' ' << 0
3242 if (n_data_sets != 0)
3244 out << n_data_sets <<
" ";
3245 for (
unsigned int i = 0; i < n_data_sets; ++i)
3250 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
3251 out << data_names[data_set]
3255 write_data(patches, n_data_sets,
true, ucd_out);
3265 template <
int dim,
int spacedim>
3269 const std::vector<std::string> &data_names,
3271 std::tuple<
unsigned int,
3283 #ifndef DEAL_II_WITH_MPI
3292 if (patches.empty())
3296 DXStream dx_out(out, flags);
3299 unsigned int offset = 0;
3301 const unsigned int n_data_sets = data_names.size();
3304 unsigned int n_nodes;
3306 std::tie(n_nodes,
n_cells) = count_nodes_and_cells(patches);
3309 out <<
"object \"vertices\" class array type float rank 1 shape "
3310 << spacedim <<
" items " << n_nodes;
3314 out <<
" lsb ieee data 0" <<
'\n';
3315 offset += n_nodes * spacedim *
sizeof(float);
3319 out <<
" data follows" <<
'\n';
3328 out <<
"object \"cells\" class array type int rank 1 shape "
3333 out <<
" lsb binary data " << offset <<
'\n';
3334 offset +=
n_cells *
sizeof(int);
3338 out <<
" data follows" <<
'\n';
3344 out <<
"attribute \"element type\" string \"";
3351 out <<
"\"" <<
'\n' <<
"attribute \"ref\" string \"positions\"" <<
'\n';
3358 out <<
"object \"neighbors\" class array type int rank 1 shape "
3362 for (
const auto &patch : patches)
3365 const unsigned int n1 = (dim > 0) ? n : 1;
3366 const unsigned int n2 = (dim > 1) ? n : 1;
3367 const unsigned int n3 = (dim > 2) ? n : 1;
3368 const unsigned int x_minus = (dim > 0) ? 0 : 0;
3369 const unsigned int x_plus = (dim > 0) ? 1 : 0;
3370 const unsigned int y_minus = (dim > 1) ? 2 : 0;
3371 const unsigned int y_plus = (dim > 1) ? 3 : 0;
3372 const unsigned int z_minus = (dim > 2) ? 4 : 0;
3373 const unsigned int z_plus = (dim > 2) ? 5 : 0;
3374 unsigned int cells_per_patch = Utilities::fixed_power<dim>(n);
3375 unsigned int dx = 1;
3376 unsigned int dy = n;
3377 unsigned int dz = n * n;
3379 const unsigned int patch_start =
3382 for (
unsigned int i3 = 0; i3 < n3; ++i3)
3383 for (
unsigned int i2 = 0; i2 < n2; ++i2)
3384 for (
unsigned int i1 = 0; i1 < n1; ++i1)
3386 const unsigned int nx = i1 *
dx;
3387 const unsigned int ny = i2 * dy;
3388 const unsigned int nz = i3 * dz;
3400 const unsigned int nn = patch.
neighbors[x_minus];
3404 << (nn * cells_per_patch + ny + nz +
dx * (n - 1));
3410 out <<
'\t' << patch_start + nx -
dx + ny + nz;
3415 const unsigned int nn = patch.
neighbors[x_plus];
3418 out << (nn * cells_per_patch + ny + nz);
3424 out <<
'\t' << patch_start + nx +
dx + ny + nz;
3431 const unsigned int nn = patch.
neighbors[y_minus];
3435 << (nn * cells_per_patch + nx + nz + dy * (n - 1));
3441 out <<
'\t' << patch_start + nx + ny - dy + nz;
3446 const unsigned int nn = patch.
neighbors[y_plus];
3449 out << (nn * cells_per_patch + nx + nz);
3455 out <<
'\t' << patch_start + nx + ny + dy + nz;
3463 const unsigned int nn = patch.
neighbors[z_minus];
3467 << (nn * cells_per_patch + nx + ny + dz * (n - 1));
3473 out <<
'\t' << patch_start + nx + ny + nz - dz;
3478 const unsigned int nn = patch.
neighbors[z_plus];
3481 out << (nn * cells_per_patch + nx + ny);
3487 out <<
'\t' << patch_start + nx + ny + nz + dz;
3495 if (n_data_sets != 0)
3497 out <<
"object \"data\" class array type float rank 1 shape "
3498 << n_data_sets <<
" items " << n_nodes;
3502 out <<
" lsb ieee data " << offset <<
'\n';
3503 offset += n_data_sets * n_nodes *
3504 ((flags.
data_double) ?
sizeof(
double) :
sizeof(float));
3508 out <<
" data follows" <<
'\n';
3513 out <<
"attribute \"dep\" string \"positions\"" <<
'\n';
3517 out <<
"object \"data\" class constantarray type float rank 0 items "
3518 << n_nodes <<
" data follows" <<
'\n'
3524 out <<
"object \"deal data\" class field" <<
'\n'
3525 <<
"component \"positions\" value \"vertices\"" <<
'\n'
3526 <<
"component \"connections\" value \"cells\"" <<
'\n'
3527 <<
"component \"data\" value \"data\"" <<
'\n';
3530 out <<
"component \"neighbors\" value \"neighbors\"" <<
'\n';
3537 out <<
"end" <<
'\n';
3555 template <
int dim,
int spacedim>
3559 const std::vector<std::string> &data_names,
3561 std::tuple<
unsigned int,
3570 #ifndef DEAL_II_WITH_MPI
3580 if (patches.empty())
3584 const unsigned int n_data_sets = data_names.size();
3588 out <<
"# This file was generated by the deal.II library." <<
'\n'
3592 <<
"# For a description of the GNUPLOT format see the GNUPLOT manual."
3599 for (
unsigned int spacedim_n = 0; spacedim_n < spacedim; ++spacedim_n)
3604 for (
const auto &data_name : data_names)
3605 out <<
'<' << data_name <<
"> ";
3611 for (
const auto &patch : patches)
3614 const unsigned int n_points_per_direction = n_subdivisions + 1;
3616 Assert((patch.
data.n_rows() == n_data_sets &&
3618 (patch.
data.n_rows() == n_data_sets + spacedim &&
3621 (n_data_sets + spacedim) :
3623 patch.
data.n_rows()));
3625 auto output_point_data =
3626 [&out, &patch, n_data_sets](
const unsigned int point_index)
mutable {
3627 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
3628 out << patch.data(data_set, point_index) <<
' ';
3637 Assert(patch.data.n_cols() == 1,
3639 n_subdivisions + 1));
3643 out << get_equispaced_location(patch, {}, n_subdivisions)
3645 output_point_data(0);
3655 Assert(patch.data.n_cols() ==
3656 Utilities::fixed_power<dim>(n_points_per_direction),
3658 n_subdivisions + 1));
3660 for (
unsigned int i1 = 0; i1 < n_points_per_direction; ++i1)
3663 out << get_equispaced_location(patch, {i1}, n_subdivisions)
3666 output_point_data(i1);
3679 Assert(patch.data.n_cols() == Utilities::fixed_power<dim>(
3680 n_points_per_direction),
3682 n_subdivisions + 1));
3684 for (
unsigned int i2 = 0; i2 < n_points_per_direction; ++i2)
3686 for (
unsigned int i1 = 0; i1 < n_points_per_direction;
3690 out << get_equispaced_location(patch,
3695 output_point_data(i1 + i2 * n_points_per_direction);
3720 out << get_node_location(patch, 0) <<
' ';
3721 output_point_data(0);
3724 out << get_node_location(patch, 1) <<
' ';
3725 output_point_data(1);
3729 out << get_node_location(patch, 2) <<
' ';
3730 output_point_data(2);
3733 out << get_node_location(patch, 2) <<
' ';
3734 output_point_data(2);
3754 Assert(patch.data.n_cols() == Utilities::fixed_power<dim>(
3755 n_points_per_direction),
3757 n_subdivisions + 1));
3762 for (
unsigned int i3 = 0; i3 < n_points_per_direction; ++i3)
3763 for (
unsigned int i2 = 0; i2 < n_points_per_direction;
3765 for (
unsigned int i1 = 0; i1 < n_points_per_direction;
3770 get_equispaced_location(patch,
3774 if (i1 < n_subdivisions)
3777 out << this_point <<
' ';
3778 output_point_data(i1 +
3779 i2 * n_points_per_direction +
3780 i3 * n_points_per_direction *
3781 n_points_per_direction);
3785 out << get_equispaced_location(patch,
3790 output_point_data((i1 + 1) +
3791 i2 * n_points_per_direction +
3792 i3 * n_points_per_direction *
3793 n_points_per_direction);
3797 out <<
'\n' <<
'\n';
3801 if (i2 < n_subdivisions)
3804 out << this_point <<
' ';
3805 output_point_data(i1 +
3806 i2 * n_points_per_direction +
3807 i3 * n_points_per_direction *
3808 n_points_per_direction);
3812 out << get_equispaced_location(patch,
3818 i1 + (i2 + 1) * n_points_per_direction +
3819 i3 * n_points_per_direction *
3820 n_points_per_direction);
3824 out <<
'\n' <<
'\n';
3828 if (i3 < n_subdivisions)
3831 out << this_point <<
' ';
3832 output_point_data(i1 +
3833 i2 * n_points_per_direction +
3834 i3 * n_points_per_direction *
3835 n_points_per_direction);
3839 out << get_equispaced_location(patch,
3845 i1 + i2 * n_points_per_direction +
3846 (i3 + 1) * n_points_per_direction *
3847 n_points_per_direction);
3850 out <<
'\n' <<
'\n';
3859 for (
const unsigned int v : {0, 1, 2, 0, 3, 2})
3861 out << get_node_location(patch, v) <<
' ';
3862 output_point_data(v);
3867 for (
const unsigned int v : {3, 1})
3869 out << get_node_location(patch, v) <<
' ';
3870 output_point_data(v);
3880 for (
const unsigned int v : {0, 1, 3, 2, 0, 4, 1})
3882 out << get_node_location(patch, v) <<
' ';
3883 output_point_data(v);
3888 for (
const unsigned int v : {2, 4, 3})
3890 out << get_node_location(patch, v) <<
' ';
3891 output_point_data(v);
3905 for (
const unsigned int v : {0, 1, 2, 0, 3, 4, 5, 3})
3907 out << get_node_location(patch, v) <<
' ';
3908 output_point_data(v);
3913 for (
const unsigned int v : {1, 4})
3915 out << get_node_location(patch, v) <<
' ';
3916 output_point_data(v);
3921 for (
const unsigned int v : {2, 5})
3923 out << get_node_location(patch, v) <<
' ';
3924 output_point_data(v);
3949 template <
int dim,
int spacedim>
3951 do_write_povray(
const std::vector<Patch<dim, spacedim>> &,
3952 const std::vector<std::string> &,
3953 const PovrayFlags &,
3957 ExcMessage(
"Writing files in POVRAY format is only supported "
3958 "for two-dimensional meshes."));
3964 do_write_povray(
const std::vector<Patch<2, 2>> &patches,
3965 const std::vector<std::string> &data_names,
3966 const PovrayFlags &flags,
3971 #ifndef DEAL_II_WITH_MPI
3980 if (patches.empty())
3983 constexpr
int dim = 2;
3985 constexpr
int spacedim = 2;
3987 const unsigned int n_data_sets = data_names.size();
3993 <<
"/* This file was generated by the deal.II library." <<
'\n'
3997 <<
" For a description of the POVRAY format see the POVRAY manual."
4002 out <<
"#include \"colors.inc\" " <<
'\n'
4003 <<
"#include \"textures.inc\" " <<
'\n';
4007 if (flags.external_data)
4008 out <<
"#include \"data.inc\" " <<
'\n';
4014 <<
"camera {" <<
'\n'
4015 <<
" location <1,4,-7>" <<
'\n'
4016 <<
" look_at <0,0,0>" <<
'\n'
4017 <<
" angle 30" <<
'\n'
4022 <<
"light_source {" <<
'\n'
4023 <<
" <1,4,-7>" <<
'\n'
4024 <<
" color Grey" <<
'\n'
4027 <<
"light_source {" <<
'\n'
4028 <<
" <0,20,0>" <<
'\n'
4029 <<
" color White" <<
'\n'
4036 double hmin = patches[0].data(0, 0);
4037 double hmax = patches[0].data(0, 0);
4039 for (
const auto &patch : patches)
4043 Assert((patch.
data.n_rows() == n_data_sets &&
4045 (patch.
data.n_rows() == n_data_sets + spacedim &&
4048 (n_data_sets + spacedim) :
4050 patch.
data.n_rows()));
4052 Utilities::fixed_power<dim>(n_subdivisions + 1),
4054 n_subdivisions + 1));
4056 for (
unsigned int i = 0; i < n_subdivisions + 1; ++i)
4057 for (
unsigned int j = 0; j < n_subdivisions + 1; ++j)
4059 const int dl = i * (n_subdivisions + 1) + j;
4060 if (patch.
data(0, dl) < hmin)
4061 hmin = patch.
data(0, dl);
4062 if (patch.
data(0, dl) > hmax)
4063 hmax = patch.
data(0, dl);
4067 out <<
"#declare HMIN=" << hmin <<
";" <<
'\n'
4068 <<
"#declare HMAX=" << hmax <<
";" <<
'\n'
4071 if (!flags.external_data)
4074 out <<
"#declare Tex=texture{" <<
'\n'
4075 <<
" pigment {" <<
'\n'
4076 <<
" gradient y" <<
'\n'
4077 <<
" scale y*(HMAX-HMIN)*" << 0.1 <<
'\n'
4078 <<
" color_map {" <<
'\n'
4079 <<
" [0.00 color Light_Purple] " <<
'\n'
4080 <<
" [0.95 color Light_Purple] " <<
'\n'
4081 <<
" [1.00 color White] " <<
'\n'
4086 if (!flags.bicubic_patch)
4089 out <<
'\n' <<
"mesh {" <<
'\n';
4093 for (
const auto &patch : patches)
4096 const unsigned int n = n_subdivisions + 1;
4097 const unsigned int d1 = 1;
4098 const unsigned int d2 = n;
4100 Assert((patch.
data.n_rows() == n_data_sets &&
4102 (patch.
data.n_rows() == n_data_sets + spacedim &&
4105 (n_data_sets + spacedim) :
4107 patch.
data.n_rows()));
4108 Assert(patch.
data.n_cols() == Utilities::fixed_power<dim>(n),
4110 n_subdivisions + 1));
4113 std::vector<Point<spacedim>> ver(n * n);
4115 for (
unsigned int i2 = 0; i2 < n; ++i2)
4116 for (
unsigned int i1 = 0; i1 < n; ++i1)
4119 ver[i1 * d1 + i2 * d2] =
4120 get_equispaced_location(patch, {i1, i2}, n_subdivisions);
4124 if (!flags.bicubic_patch)
4127 std::vector<Point<3>> nrml;
4137 for (
unsigned int i = 0; i < n; ++i)
4138 for (
unsigned int j = 0; j < n; ++j)
4140 const unsigned int il = (i == 0) ? i : (i - 1);
4141 const unsigned int ir =
4142 (i == n_subdivisions) ? i : (i + 1);
4143 const unsigned int jl = (j == 0) ? j : (j - 1);
4144 const unsigned int jr =
4145 (j == n_subdivisions) ? j : (j + 1);
4148 ver[ir * d1 + j * d2](0) - ver[il * d1 + j * d2](0);
4149 h1(1) = patch.
data(0, ir * d1 + j * d2) -
4150 patch.
data(0, il * d1 + j * d2);
4152 ver[ir * d1 + j * d2](1) - ver[il * d1 + j * d2](1);
4155 ver[i * d1 + jr * d2](0) - ver[i * d1 + jl * d2](0);
4156 h2(1) = patch.
data(0, i * d1 + jr * d2) -
4157 patch.
data(0, i * d1 + jl * d2);
4159 ver[i * d1 + jr * d2](1) - ver[i * d1 + jl * d2](1);
4161 nrml[i * d1 + j * d2](0) =
4162 h1(1) * h2(2) - h1(2) * h2(1);
4163 nrml[i * d1 + j * d2](1) =
4164 h1(2) * h2(0) - h1(0) * h2(2);
4165 nrml[i * d1 + j * d2](2) =
4166 h1(0) * h2(1) - h1(1) * h2(0);
4170 std::sqrt(std::pow(nrml[i * d1 + j * d2](0), 2.) +
4171 std::pow(nrml[i * d1 + j * d2](1), 2.) +
4172 std::pow(nrml[i * d1 + j * d2](2), 2.));
4174 if (nrml[i * d1 + j * d2](1) < 0)
4177 for (
unsigned int k = 0; k < 3; ++k)
4178 nrml[i * d1 + j * d2](k) /=
norm;
4183 for (
unsigned int i = 0; i < n_subdivisions; ++i)
4184 for (
unsigned int j = 0; j < n_subdivisions; ++j)
4187 const int dl = i * d1 + j * d2;
4193 out <<
"smooth_triangle {" <<
'\n'
4194 <<
"\t<" << ver[dl](0) <<
"," << patch.
data(0, dl)
4195 <<
"," << ver[dl](1) <<
">, <" << nrml[dl](0)
4196 <<
", " << nrml[dl](1) <<
", " << nrml[dl](2)
4198 out <<
" \t<" << ver[dl + d1](0) <<
","
4199 << patch.
data(0, dl + d1) <<
"," << ver[dl + d1](1)
4200 <<
">, <" << nrml[dl + d1](0) <<
", "
4201 << nrml[dl + d1](1) <<
", " << nrml[dl + d1](2)
4203 out <<
"\t<" << ver[dl + d1 + d2](0) <<
","
4204 << patch.
data(0, dl + d1 + d2) <<
","
4205 << ver[dl + d1 + d2](1) <<
">, <"
4206 << nrml[dl + d1 + d2](0) <<
", "
4207 << nrml[dl + d1 + d2](1) <<
", "
4208 << nrml[dl + d1 + d2](2) <<
">}" <<
'\n';
4211 out <<
"smooth_triangle {" <<
'\n'
4212 <<
"\t<" << ver[dl](0) <<
"," << patch.
data(0, dl)
4213 <<
"," << ver[dl](1) <<
">, <" << nrml[dl](0)
4214 <<
", " << nrml[dl](1) <<
", " << nrml[dl](2)
4216 out <<
"\t<" << ver[dl + d1 + d2](0) <<
","
4217 << patch.
data(0, dl + d1 + d2) <<
","
4218 << ver[dl + d1 + d2](1) <<
">, <"
4219 << nrml[dl + d1 + d2](0) <<
", "
4220 << nrml[dl + d1 + d2](1) <<
", "
4221 << nrml[dl + d1 + d2](2) <<
">," <<
'\n';
4222 out <<
"\t<" << ver[dl + d2](0) <<
","
4223 << patch.
data(0, dl + d2) <<
"," << ver[dl + d2](1)
4224 <<
">, <" << nrml[dl + d2](0) <<
", "
4225 << nrml[dl + d2](1) <<
", " << nrml[dl + d2](2)
4231 out <<
"triangle {" <<
'\n'
4232 <<
"\t<" << ver[dl](0) <<
"," << patch.
data(0, dl)
4233 <<
"," << ver[dl](1) <<
">," <<
'\n';
4234 out <<
"\t<" << ver[dl + d1](0) <<
","
4235 << patch.
data(0, dl + d1) <<
"," << ver[dl + d1](1)
4237 out <<
"\t<" << ver[dl + d1 + d2](0) <<
","
4238 << patch.
data(0, dl + d1 + d2) <<
","
4239 << ver[dl + d1 + d2](1) <<
">}" <<
'\n';
4242 out <<
"triangle {" <<
'\n'
4243 <<
"\t<" << ver[dl](0) <<
"," << patch.
data(0, dl)
4244 <<
"," << ver[dl](1) <<
">," <<
'\n';
4245 out <<
"\t<" << ver[dl + d1 + d2](0) <<
","
4246 << patch.
data(0, dl + d1 + d2) <<
","
4247 << ver[dl + d1 + d2](1) <<
">," <<
'\n';
4248 out <<
"\t<" << ver[dl + d2](0) <<
","
4249 << patch.
data(0, dl + d2) <<
"," << ver[dl + d2](1)
4257 Assert(n_subdivisions == 3,
4260 <<
"bicubic_patch {" <<
'\n'
4261 <<
" type 0" <<
'\n'
4262 <<
" flatness 0" <<
'\n'
4263 <<
" u_steps 0" <<
'\n'
4264 <<
" v_steps 0" <<
'\n';
4265 for (
int i = 0; i < 16; ++i)
4267 out <<
"\t<" << ver[i](0) <<
"," << patch.
data(0, i) <<
","
4268 << ver[i](1) <<
">";
4273 out <<
" texture {Tex}" <<
'\n' <<
"}" <<
'\n';
4277 if (!flags.bicubic_patch)
4280 out <<
" texture {Tex}" <<
'\n' <<
"}" <<
'\n' <<
'\n';
4292 template <
int dim,
int spacedim>
4296 const std::vector<std::string> &data_names,
4298 std::tuple<
unsigned int,
4305 do_write_povray(patches, data_names, flags, out);
4310 template <
int dim,
int spacedim>
4314 const std::vector<std::string> & ,
4316 std::tuple<
unsigned int,
4328 template <
int spacedim>
4332 const std::vector<std::string> & ,
4334 std::tuple<
unsigned int,
4343 #ifndef DEAL_II_WITH_MPI
4352 if (patches.empty())
4362 std::multiset<EpsCell2d> cells;
4371 double heights[4] = {0, 0, 0, 0};
4375 for (
const auto &patch : patches)
4378 const unsigned int n = n_subdivisions + 1;
4379 const unsigned int d1 = 1;
4380 const unsigned int d2 = n;
4382 for (
unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
4383 for (
unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
4387 get_equispaced_location(patch, {i1, i2}, n_subdivisions);
4389 get_equispaced_location(patch, {i1 + 1, i2}, n_subdivisions);
4391 get_equispaced_location(patch, {i1, i2 + 1}, n_subdivisions);
4392 points[3] = get_equispaced_location(patch,
4400 patch.
data.n_rows() == 0,
4403 patch.
data.n_rows()));
4405 patch.
data.n_rows() != 0 ?
4409 heights[1] = patch.
data.n_rows() != 0 ?
4411 (i1 + 1) * d1 + i2 * d2) *
4414 heights[2] = patch.
data.n_rows() != 0 ?
4416 i1 * d1 + (i2 + 1) * d2) *
4419 heights[3] = patch.
data.n_rows() != 0 ?
4421 (i1 + 1) * d1 + (i2 + 1) * d2) *
4428 for (
unsigned int i = 0; i < 4; ++i)
4429 heights[i] = points[i](2);
4450 cz = -std::cos(flags.
turn_angle * 2 * pi / 360.),
4453 sz = std::sin(flags.
turn_angle * 2 * pi / 360.);
4454 for (
unsigned int vertex = 0; vertex < 4; ++vertex)
4456 const double x = points[vertex](0), y = points[vertex](1),
4457 z = -heights[vertex];
4459 eps_cell.vertices[vertex](0) = -cz * x + sz * y;
4460 eps_cell.vertices[vertex](1) =
4461 -cx * sz * x - cx * cz * y - sx * z;
4485 (points[0] + points[1] + points[2] + points[3]) / 4;
4486 const double center_height =
4487 -(heights[0] + heights[1] + heights[2] + heights[3]) / 4;
4490 eps_cell.depth = -sx * sz * center_point(0) -
4491 sx * cz * center_point(1) + cx * center_height;
4496 patch.
data.n_rows() == 0,
4499 patch.
data.n_rows()));
4500 const double color_values[4] = {
4501 patch.
data.n_rows() != 0 ?
4505 patch.
data.n_rows() != 0 ?
4509 patch.
data.n_rows() != 0 ?
4513 patch.
data.n_rows() != 0 ?
4515 (i1 + 1) * d1 + (i2 + 1) * d2) :
4519 eps_cell.color_value = (color_values[0] + color_values[1] +
4520 color_values[3] + color_values[2]) /
4525 std::min(min_color_value, eps_cell.color_value);
4527 std::max(max_color_value, eps_cell.color_value);
4531 cells.insert(eps_cell);
4537 double x_min = cells.begin()->vertices[0](0);
4538 double x_max = x_min;
4539 double y_min = cells.begin()->vertices[0](1);
4540 double y_max = y_min;
4542 for (
const auto &cell : cells)
4543 for (
const auto &vertex : cell.vertices)
4545 x_min =
std::min(x_min, vertex(0));
4546 x_max =
std::max(x_max, vertex(0));
4547 y_min =
std::min(y_min, vertex(1));
4548 y_max =
std::max(y_max, vertex(1));
4553 const double scale =
4557 const Point<2> offset(x_min, y_min);
4562 out <<
"%!PS-Adobe-2.0 EPSF-1.2" <<
'\n'
4563 <<
"%%Title: deal.II Output" <<
'\n'
4564 <<
"%%Creator: the deal.II library" <<
'\n'
4567 <<
"%%BoundingBox: "
4571 <<
static_cast<unsigned int>((x_max - x_min) *
scale + 0.5) <<
' '
4572 <<
static_cast<unsigned int>((y_max - y_min) *
scale + 0.5) <<
'\n';
4581 out <<
"/m {moveto} bind def" <<
'\n'
4582 <<
"/l {lineto} bind def" <<
'\n'
4583 <<
"/s {setrgbcolor} bind def" <<
'\n'
4584 <<
"/sg {setgray} bind def" <<
'\n'
4585 <<
"/lx {lineto closepath stroke} bind def" <<
'\n'
4586 <<
"/lf {lineto closepath fill} bind def" <<
'\n';
4588 out <<
"%%EndProlog" <<
'\n' <<
'\n';
4590 out << flags.
line_width <<
" setlinewidth" <<
'\n';
4598 if (max_color_value == min_color_value)
4599 max_color_value = min_color_value + 1;
4603 for (
const auto &cell : cells)
4616 out << rgb_values.
red <<
" sg ";
4618 out << rgb_values.
red <<
' ' << rgb_values.
green <<
' '
4619 << rgb_values.
blue <<
" s ";
4624 out << (cell.vertices[0] - offset) *
scale <<
" m "
4625 << (cell.vertices[1] - offset) *
scale <<
" l "
4626 << (cell.vertices[3] - offset) *
scale <<
" l "
4627 << (cell.vertices[2] - offset) *
scale <<
" lf" <<
'\n';
4632 << (cell.vertices[0] - offset) *
scale <<
" m "
4633 << (cell.vertices[1] - offset) *
scale <<
" l "
4634 << (cell.vertices[3] - offset) *
scale <<
" l "
4635 << (cell.vertices[2] - offset) *
scale <<
" lx" <<
'\n';
4637 out <<
"showpage" <<
'\n';
4646 template <
int dim,
int spacedim>
4650 const std::vector<std::string> &data_names,
4652 std::tuple<
unsigned int,
4668 #ifndef DEAL_II_WITH_MPI
4677 if (patches.empty())
4681 GmvStream gmv_out(out, flags);
4682 const unsigned int n_data_sets = data_names.size();
4685 Assert((patches[0].data.n_rows() == n_data_sets &&
4686 !patches[0].points_are_available) ||
4687 (patches[0].data.n_rows() == n_data_sets + spacedim &&
4688 patches[0].points_are_available),
4690 (n_data_sets + spacedim) :
4692 patches[0].data.n_rows()));
4696 out <<
"gmvinput ascii" <<
'\n' <<
'\n';
4699 unsigned int n_nodes;
4701 std::tie(n_nodes,
n_cells) = count_nodes_and_cells(patches);
4716 [&patches]() {
return create_global_data_table(patches); });
4722 out <<
"nodes " << n_nodes <<
'\n';
4723 for (
unsigned int d = 0;
d < spacedim; ++
d)
4725 gmv_out.selected_component =
d;
4731 for (
unsigned int d = spacedim;
d < 3; ++
d)
4733 for (
unsigned int i = 0; i < n_nodes; ++i)
4740 out <<
"cells " <<
n_cells <<
'\n';
4745 out <<
"variable" <<
'\n';
4749 std::move(*create_global_data_table_task.
return_value());
4753 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
4755 out << data_names[data_set] <<
" 1" <<
'\n';
4757 data_vectors[data_set].
end(),
4758 std::ostream_iterator<double>(out,
" "));
4759 out <<
'\n' <<
'\n';
4765 out <<
"endvars" <<
'\n';
4768 out <<
"endgmv" <<
'\n';
4779 template <
int dim,
int spacedim>
4783 const std::vector<std::string> &data_names,
4785 std::tuple<
unsigned int,
4799 #ifndef DEAL_II_WITH_MPI
4808 if (patches.empty())
4812 TecplotStream tecplot_out(out, flags);
4814 const unsigned int n_data_sets = data_names.size();
4817 Assert((patches[0].data.n_rows() == n_data_sets &&
4818 !patches[0].points_are_available) ||
4819 (patches[0].data.n_rows() == n_data_sets + spacedim &&
4820 patches[0].points_are_available),
4822 (n_data_sets + spacedim) :
4824 patches[0].data.n_rows()));
4827 unsigned int n_nodes;
4829 std::tie(n_nodes,
n_cells) = count_nodes_and_cells(patches);
4835 <<
"# This file was generated by the deal.II library." <<
'\n'
4839 <<
"# For a description of the Tecplot format see the Tecplot documentation."
4844 out <<
"Variables=";
4852 out <<
"\"x\", \"y\"";
4855 out <<
"\"x\", \"y\", \"z\"";
4861 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
4862 out <<
", \"" << data_names[data_set] <<
"\"";
4868 out <<
"t=\"" << flags.
zone_name <<
"\" ";
4871 out <<
"strandid=1, solutiontime=" << flags.
solution_time <<
", ";
4873 out <<
"f=feblock, n=" << n_nodes <<
", e=" <<
n_cells
4874 <<
", et=" << tecplot_cell_type[dim] <<
'\n';
4891 [&patches]() {
return create_global_data_table(patches); });
4897 for (
unsigned int d = 0;
d < spacedim; ++
d)
4899 tecplot_out.selected_component =
d;
4910 std::move(*create_global_data_table_task.
return_value());
4913 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
4916 data_vectors[data_set].
end(),
4917 std::ostream_iterator<double>(out,
"\n"));
4932 template <
int dim,
int spacedim>
4936 const std::vector<std::string> &data_names,
4938 std::tuple<
unsigned int,
4942 &nonscalar_data_ranges,
4948 #ifndef DEAL_II_WITH_MPI
4957 if (patches.empty())
4961 VtkStream vtk_out(out, flags);
4963 const unsigned int n_data_sets = data_names.size();
4965 if (patches[0].points_are_available)
4977 out <<
"# vtk DataFile Version 3.0" <<
'\n'
4978 <<
"#This file was generated by the deal.II library";
4986 out <<
'\n' <<
"ASCII" <<
'\n';
4988 out <<
"DATASET UNSTRUCTURED_GRID\n" <<
'\n';
4995 const unsigned int n_metadata =
5000 out <<
"FIELD FieldData " << n_metadata <<
'\n';
5004 out <<
"CYCLE 1 1 int\n" << flags.
cycle <<
'\n';
5008 out <<
"TIME 1 1 double\n" << flags.
time <<
'\n';
5014 unsigned int n_nodes;
5016 unsigned int n_points_and_n_cells;
5017 std::tie(n_nodes,
n_cells, n_points_and_n_cells) =
5033 [&patches]() {
return create_global_data_table(patches); });
5039 out <<
"POINTS " << n_nodes <<
" double" <<
'\n';
5044 out <<
"CELLS " <<
n_cells <<
' ' << n_points_and_n_cells <<
'\n';
5052 out <<
"CELL_TYPES " <<
n_cells <<
'\n';
5056 for (
const auto &patch : patches)
5058 const auto vtk_cell_id =
5061 for (
unsigned int i = 0; i < vtk_cell_id[1]; ++i)
5062 out <<
' ' << vtk_cell_id[0];
5071 std::move(*create_global_data_table_task.
return_value());
5076 out <<
"POINT_DATA " << n_nodes <<
'\n';
5080 std::vector<bool> data_set_written(n_data_sets,
false);
5081 for (
const auto &nonscalar_data_range : nonscalar_data_ranges)
5086 "The VTK writer does not currently support outputting "
5087 "tensor data. Use the VTU writer instead."));
5090 std::get<0>(nonscalar_data_range),
5092 std::get<0>(nonscalar_data_range)));
5093 AssertThrow(std::get<1>(nonscalar_data_range) < n_data_sets,
5097 AssertThrow(std::get<1>(nonscalar_data_range) + 1 -
5098 std::get<0>(nonscalar_data_range) <=
5101 "Can't declare a vector with more than 3 components "
5105 for (
unsigned int i = std::get<0>(nonscalar_data_range);
5106 i <= std::get<1>(nonscalar_data_range);
5108 data_set_written[i] =
true;
5114 if (!std::get<2>(nonscalar_data_range).empty())
5115 out << std::get<2>(nonscalar_data_range);
5118 for (
unsigned int i = std::get<0>(nonscalar_data_range);
5119 i < std::get<1>(nonscalar_data_range);
5121 out << data_names[i] <<
"__";
5122 out << data_names[std::get<1>(nonscalar_data_range)];
5125 out <<
" double" <<
'\n';
5128 for (
unsigned int n = 0; n < n_nodes; ++n)
5130 switch (std::get<1>(nonscalar_data_range) -
5131 std::get<0>(nonscalar_data_range))
5134 out << data_vectors(std::get<0>(nonscalar_data_range), n)
5139 out << data_vectors(std::get<0>(nonscalar_data_range), n)
5141 << data_vectors(std::get<0>(nonscalar_data_range) + 1, n)
5145 out << data_vectors(std::get<0>(nonscalar_data_range), n)
5147 << data_vectors(std::get<0>(nonscalar_data_range) + 1, n)
5149 << data_vectors(std::get<0>(nonscalar_data_range) + 2, n)
5162 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
5163 if (data_set_written[data_set] ==
false)
5165 out <<
"SCALARS " << data_names[data_set] <<
" double 1" <<
'\n'
5166 <<
"LOOKUP_TABLE default" <<
'\n';
5168 data_vectors[data_set].
end(),
5169 std::ostream_iterator<double>(out,
" "));
5185 out <<
"<?xml version=\"1.0\" ?> \n";
5187 out <<
"# vtk DataFile Version 3.0" <<
'\n'
5188 <<
"#This file was generated by the deal.II library";
5199 out <<
"<VTKFile type=\"UnstructuredGrid\" version=\"2.2\"";
5201 out <<
"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\"";
5202 if (deal_ii_with_zlib &&
5204 out <<
" compressor=\"vtkZLibDataCompressor\"";
5205 #ifdef DEAL_II_WORDS_BIGENDIAN
5206 out <<
" byte_order=\"BigEndian\"";
5208 out <<
" byte_order=\"LittleEndian\"";
5212 out <<
"<UnstructuredGrid>";
5222 out <<
" </UnstructuredGrid>\n";
5223 out <<
"</VTKFile>\n";
5228 template <
int dim,
int spacedim>
5232 const std::vector<std::string> &data_names,
5234 std::tuple<
unsigned int,
5238 &nonscalar_data_ranges,
5243 write_vtu_main(patches, data_names, nonscalar_data_ranges, flags, out);
5250 template <
int dim,
int spacedim>
5254 const std::vector<std::string> &data_names,
5256 std::tuple<
unsigned int,
5260 &nonscalar_data_ranges,
5274 unit.second.find(
'\"') == std::string::npos,
5276 "A physical unit you provided, <" + unit.second +
5277 ">, contained a quotation mark character. This is not allowed."));
5280 #ifndef DEAL_II_WITH_MPI
5289 if (patches.empty())
5294 out <<
"<Piece NumberOfPoints=\"0\" NumberOfCells=\"0\" >\n"
5296 <<
"<DataArray type=\"UInt8\" Name=\"types\"></DataArray>\n"
5298 <<
" <PointData Scalars=\"scalars\">\n";
5299 std::vector<bool> data_set_written(data_names.size(),
false);
5300 for (
const auto &nonscalar_data_range : nonscalar_data_ranges)
5303 for (
unsigned int i = std::get<0>(nonscalar_data_range);
5304 i <= std::get<1>(nonscalar_data_range);
5306 data_set_written[i] =
true;
5310 out <<
" <DataArray type=\"Float32\" Name=\"";
5312 if (!std::get<2>(nonscalar_data_range).empty())
5313 out << std::get<2>(nonscalar_data_range);
5316 for (
unsigned int i = std::get<0>(nonscalar_data_range);
5317 i < std::get<1>(nonscalar_data_range);
5319 out << data_names[i] <<
"__";
5320 out << data_names[std::get<1>(nonscalar_data_range)];
5323 out <<
"\" NumberOfComponents=\"3\"></DataArray>\n";
5326 for (
unsigned int data_set = 0; data_set < data_names.size();
5328 if (data_set_written[data_set] ==
false)
5330 out <<
" <DataArray type=\"Float32\" Name=\""
5331 << data_names[data_set] <<
"\"></DataArray>\n";
5334 out <<
" </PointData>\n";
5335 out <<
"</Piece>\n";
5349 const unsigned int n_metadata =
5353 out <<
"<FieldData>\n";
5358 <<
"<DataArray type=\"Float32\" Name=\"CYCLE\" NumberOfTuples=\"1\" format=\"ascii\">"
5359 << flags.
cycle <<
"</DataArray>\n";
5364 <<
"<DataArray type=\"Float32\" Name=\"TIME\" NumberOfTuples=\"1\" format=\"ascii\">"
5365 << flags.
time <<
"</DataArray>\n";
5369 out <<
"</FieldData>\n";
5373 const unsigned int n_data_sets = data_names.size();
5376 if (patches[0].points_are_available)
5385 const char *ascii_or_binary =
5386 (deal_ii_with_zlib &&
5393 unsigned int n_nodes;
5395 std::tie(n_nodes,
n_cells, std::ignore) =
5403 const auto stringize_vertex_information = [&patches,
5407 ascii_or_binary]() {
5408 std::ostringstream o;
5410 o <<
" <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\""
5411 << ascii_or_binary <<
"\">\n";
5412 const std::vector<Point<spacedim>> node_positions =
5417 std::vector<float> node_coordinates_3d;
5418 node_coordinates_3d.reserve(node_positions.size() * 3);
5419 for (
const auto &node_position : node_positions)
5421 for (
unsigned int d = 0;
d < 3; ++
d)
5423 node_coordinates_3d.emplace_back(node_position[
d]);
5425 node_coordinates_3d.emplace_back(0.0f);
5427 o << vtu_stringize_array(node_coordinates_3d,
5431 o <<
" </DataArray>\n";
5432 o <<
" </Points>\n\n";
5441 const auto stringize_cell_to_vertex_information = [&patches,
5445 out.precision()]() {
5446 std::ostringstream o;
5449 o <<
" <DataArray type=\"Int32\" Name=\"connectivity\" format=\""
5450 << ascii_or_binary <<
"\">\n";
5452 std::vector<int32_t> cells;
5455 unsigned int first_vertex_of_patch = 0;
5457 for (
const auto &patch : patches)
5469 const unsigned int n_points = patch.
data.n_cols();
5470 Assert((dim == 2 && n_points == 6) ||
5471 (dim == 3 && n_points == 10),
5474 if (deal_ii_with_zlib &&
5478 for (
unsigned int i = 0; i < n_points; ++i)
5479 cells.push_back(first_vertex_of_patch + i);
5483 for (
unsigned int i = 0; i < n_points; ++i)
5484 o <<
'\t' << first_vertex_of_patch + i;
5488 first_vertex_of_patch += n_points;
5493 else if (patch.
reference_cell != ReferenceCells::get_hypercube<dim>())
5497 const unsigned int n_points = patch.
data.n_cols();
5499 if (deal_ii_with_zlib &&
5503 for (
unsigned int i = 0; i < n_points; ++i)
5505 first_vertex_of_patch +
5510 for (
unsigned int i = 0; i < n_points; ++i)
5512 << (first_vertex_of_patch +
5517 first_vertex_of_patch += n_points;
5522 const unsigned int n_points_per_direction = n_subdivisions + 1;
5524 std::vector<unsigned> local_vertex_order;
5528 const auto flush_current_cell = [&flags,
5531 first_vertex_of_patch,
5532 &local_vertex_order]() {
5533 if (deal_ii_with_zlib &&
5537 for (
const auto &c : local_vertex_order)
5538 cells.push_back(first_vertex_of_patch + c);
5542 for (
const auto &c : local_vertex_order)
5543 o <<
'\t' << first_vertex_of_patch + c;
5547 local_vertex_order.clear();
5552 local_vertex_order.reserve(Utilities::fixed_power<dim>(2));
5558 local_vertex_order.emplace_back(0);
5559 flush_current_cell();
5565 for (
unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
5567 const unsigned int starting_offset = i1;
5568 local_vertex_order.emplace_back(starting_offset);
5569 local_vertex_order.emplace_back(starting_offset +
5571 flush_current_cell();
5578 for (
unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
5579 for (
unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
5581 const unsigned int starting_offset =
5582 i2 * n_points_per_direction + i1;
5583 local_vertex_order.emplace_back(
5585 local_vertex_order.emplace_back(
5586 starting_offset + 1);
5587 local_vertex_order.emplace_back(
5588 starting_offset + n_points_per_direction + 1);
5589 local_vertex_order.emplace_back(
5590 starting_offset + n_points_per_direction);
5591 flush_current_cell();
5598 for (
unsigned int i3 = 0; i3 < n_subdivisions; ++i3)
5599 for (
unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
5600 for (
unsigned int i1 = 0; i1 < n_subdivisions;
5603 const unsigned int starting_offset =
5604 i3 * n_points_per_direction *
5605 n_points_per_direction +
5606 i2 * n_points_per_direction + i1;
5607 local_vertex_order.emplace_back(
5609 local_vertex_order.emplace_back(
5610 starting_offset + 1);
5611 local_vertex_order.emplace_back(
5612 starting_offset + n_points_per_direction +
5614 local_vertex_order.emplace_back(
5615 starting_offset + n_points_per_direction);
5616 local_vertex_order.emplace_back(
5617 starting_offset + n_points_per_direction *
5618 n_points_per_direction);
5619 local_vertex_order.emplace_back(
5621 n_points_per_direction *
5622 n_points_per_direction +
5624 local_vertex_order.emplace_back(
5626 n_points_per_direction *
5627 n_points_per_direction +
5628 n_points_per_direction + 1);
5629 local_vertex_order.emplace_back(
5631 n_points_per_direction *
5632 n_points_per_direction +
5633 n_points_per_direction);
5634 flush_current_cell();
5645 local_vertex_order.resize(
5646 Utilities::fixed_power<dim>(n_points_per_direction));
5654 "Point-like cells should not be possible "
5655 "when writing higher-order cells."));
5660 for (
unsigned int i1 = 0; i1 < n_subdivisions + 1;
5663 const unsigned int local_index = i1;
5664 const unsigned int connectivity_index =
5666 .template vtk_lexicographic_to_node_index<1>(
5670 local_vertex_order[connectivity_index] =
5672 flush_current_cell();
5679 for (
unsigned int i2 = 0; i2 < n_subdivisions + 1;
5681 for (
unsigned int i1 = 0; i1 < n_subdivisions + 1;
5684 const unsigned int local_index =
5685 i2 * n_points_per_direction + i1;
5686 const unsigned int connectivity_index =
5688 .template vtk_lexicographic_to_node_index<
5690 {{n_subdivisions, n_subdivisions}},
5692 local_vertex_order[connectivity_index] =
5695 flush_current_cell();
5701 for (
unsigned int i3 = 0; i3 < n_subdivisions + 1;
5703 for (
unsigned int i2 = 0; i2 < n_subdivisions + 1;
5705 for (
unsigned int i1 = 0; i1 < n_subdivisions + 1;
5708 const unsigned int local_index =
5709 i3 * n_points_per_direction *
5710 n_points_per_direction +
5711 i2 * n_points_per_direction + i1;
5712 const unsigned int connectivity_index =
5714 .template vtk_lexicographic_to_node_index<
5720 local_vertex_order[connectivity_index] =
5724 flush_current_cell();
5734 first_vertex_of_patch +=
5743 o << vtu_stringize_array(cells,
5748 o <<
" </DataArray>\n";
5768 const auto stringize_cell_offset_and_type_information =
5773 output_precision = out.precision()]() {
5774 std::ostringstream o;
5776 o <<
" <DataArray type=\"Int32\" Name=\"offsets\" format=\""
5777 << ascii_or_binary <<
"\">\n";
5779 std::vector<int32_t> offsets;
5784 std::vector<unsigned int> cell_types;
5787 unsigned int first_vertex_of_patch = 0;
5789 for (
const auto &patch : patches)
5791 const auto vtk_cell_id =
5794 for (
unsigned int i = 0; i < vtk_cell_id[1]; ++i)
5796 cell_types.push_back(vtk_cell_id[0]);
5797 first_vertex_of_patch += vtk_cell_id[2];
5798 offsets.push_back(first_vertex_of_patch);
5802 o << vtu_stringize_array(offsets,
5806 o <<
" </DataArray>\n";
5808 o <<
" <DataArray type=\"UInt8\" Name=\"types\" format=\""
5809 << ascii_or_binary <<
"\">\n";
5811 if (deal_ii_with_zlib &&
5814 std::vector<uint8_t> cell_types_uint8_t(cell_types.size());
5815 for (
unsigned int i = 0; i < cell_types.size(); ++i)
5816 cell_types_uint8_t[i] =
static_cast<std::uint8_t
>(cell_types[i]);
5818 o << vtu_stringize_array(cell_types_uint8_t,
5824 o << vtu_stringize_array(cell_types,
5830 o <<
" </DataArray>\n";
5840 const auto stringize_nonscalar_data_range =
5846 output_precision = out.precision()](
const Table<2, float> &data_vectors,
5847 const auto &range) {
5848 std::ostringstream o;
5850 const auto first_component = std::get<0>(range);
5851 const auto last_component = std::get<1>(range);
5852 const auto &name = std::get<2>(range);
5853 const bool is_tensor =
5854 (std::get<3>(range) ==
5856 const unsigned int n_components = (is_tensor ? 9 : 3);
5863 AssertThrow((last_component + 1 - first_component <= 9),
5865 "Can't declare a tensor with more than 9 components "
5866 "in VTK/VTU format."));
5870 AssertThrow((last_component + 1 - first_component <= 3),
5872 "Can't declare a vector with more than 3 components "
5873 "in VTK/VTU format."));
5878 o <<
" <DataArray type=\"Float32\" Name=\"";
5884 for (
unsigned int i = first_component; i < last_component; ++i)
5885 o << data_names[i] <<
"__";
5886 o << data_names[last_component];
5889 o <<
"\" NumberOfComponents=\"" << n_components <<
"\" format=\""
5890 << ascii_or_binary <<
"\"";
5909 std::vector<float> data;
5910 data.reserve(n_nodes * n_components);
5912 for (
unsigned int n = 0; n < n_nodes; ++n)
5916 switch (last_component - first_component)
5919 data.push_back(data_vectors(first_component, n));
5925 data.push_back(data_vectors(first_component, n));
5926 data.push_back(data_vectors(first_component + 1, n));
5931 data.push_back(data_vectors(first_component, n));
5932 data.push_back(data_vectors(first_component + 1, n));
5933 data.push_back(data_vectors(first_component + 2, n));
5946 const unsigned int size = last_component - first_component + 1;
5950 vtk_data[0][0] = data_vectors(first_component, n);
5955 for (
unsigned int c = 0; c < size; ++c)
5959 vtk_data[ind[0]][ind[1]] =
5960 data_vectors(first_component + c, n);
5966 for (
unsigned int c = 0; c < size; ++c)
5970 vtk_data[ind[0]][ind[1]] =
5971 data_vectors(first_component + c, n);
5982 for (
unsigned int i = 0; i < 3; ++i)
5983 for (
unsigned int j = 0; j < 3; ++j)
5984 data.push_back(vtk_data[i][j]);
5988 o << vtu_stringize_array(data,
5992 o <<
" </DataArray>\n";
5997 const auto stringize_scalar_data_set =
6001 output_precision = out.precision()](
const Table<2, float> &data_vectors,
6002 const unsigned int data_set) {
6003 std::ostringstream o;
6005 o <<
" <DataArray type=\"Float32\" Name=\"" << data_names[data_set]
6006 <<
"\" format=\"" << ascii_or_binary <<
"\"";
6010 o <<
" units=\"" << flags.
physical_units.at(data_names[data_set])
6015 const std::vector<float> data(data_vectors[data_set].
begin(),
6016 data_vectors[data_set].
end());
6017 o << vtu_stringize_array(data,
6021 o <<
" </DataArray>\n";
6040 return create_global_data_table<dim, spacedim, float>(patches);
6054 std::move(*create_global_data_table_task.
return_value());
6060 std::vector<bool> data_set_handled(n_data_sets,
false);
6061 for (
const auto &range : nonscalar_data_ranges)
6064 const auto first_component = std::get<0>(range);
6065 const auto last_component = std::get<1>(range);
6066 for (
unsigned int i = first_component; i <= last_component; ++i)
6067 data_set_handled[i] =
true;
6070 return stringize_nonscalar_data_range(data_vectors, range);
6075 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
6076 if (data_set_handled[data_set] ==
false)
6079 return stringize_scalar_data_set(data_vectors, data_set);
6085 out <<
"<Piece NumberOfPoints=\"" << n_nodes <<
"\" NumberOfCells=\""
6087 for (
const auto &s : mesh_tasks.return_values())
6089 out <<
" <PointData Scalars=\"scalars\">\n";
6092 out <<
" </PointData>\n";
6093 out <<
" </Piece>\n";
6107 const std::vector<std::string> &piece_names,
6108 const std::vector<std::string> &data_names,
6110 std::tuple<
unsigned int,
6114 &nonscalar_data_ranges,
6127 unit.second.find(
'\"') == std::string::npos,
6129 "A physical unit you provided, <" + unit.second +
6130 ">, contained a quotation mark character. This is not allowed."));
6133 const unsigned int n_data_sets = data_names.size();
6135 out <<
"<?xml version=\"1.0\"?>\n";
6138 out <<
"#This file was generated by the deal.II library"
6143 <<
"<VTKFile type=\"PUnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
6144 out <<
" <PUnstructuredGrid GhostLevel=\"0\">\n";
6145 out <<
" <PPointData Scalars=\"scalars\">\n";
6148 std::vector<bool> data_set_written(n_data_sets,
false);
6149 for (
const auto &nonscalar_data_range : nonscalar_data_ranges)
6151 const auto first_component = std::get<0>(nonscalar_data_range);
6152 const auto last_component = std::get<1>(nonscalar_data_range);
6153 const bool is_tensor =
6154 (std::get<3>(nonscalar_data_range) ==
6156 const unsigned int n_components = (is_tensor ? 9 : 3);
6163 AssertThrow((last_component + 1 - first_component <= 9),
6165 "Can't declare a tensor with more than 9 components "
6170 Assert((last_component + 1 - first_component <= 3),
6172 "Can't declare a vector with more than 3 components "
6177 for (
unsigned int i = std::get<0>(nonscalar_data_range);
6178 i <= std::get<1>(nonscalar_data_range);
6180 data_set_written[i] =
true;
6184 out <<
" <PDataArray type=\"Float32\" Name=\"";
6186 const std::string &name = std::get<2>(nonscalar_data_range);
6191 for (
unsigned int i = std::get<0>(nonscalar_data_range);
6192 i < std::get<1>(nonscalar_data_range);
6194 out << data_names[i] <<
"__";
6195 out << data_names[std::get<1>(nonscalar_data_range)];
6198 out <<
"\" NumberOfComponents=\"" << n_components
6199 <<
"\" format=\"ascii\"";
6211 data_names[std::get<1>(nonscalar_data_range)]) !=
6215 data_names[std::get<1>(nonscalar_data_range)])
6223 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
6224 if (data_set_written[data_set] ==
false)
6226 out <<
" <PDataArray type=\"Float32\" Name=\""
6227 << data_names[data_set] <<
"\" format=\"ascii\"";
6231 out <<
" units=\"" << flags.
physical_units.at(data_names[data_set])
6237 out <<
" </PPointData>\n";
6239 out <<
" <PPoints>\n";
6240 out <<
" <PDataArray type=\"Float32\" NumberOfComponents=\"3\"/>\n";
6241 out <<
" </PPoints>\n";
6243 for (
const auto &piece_name : piece_names)
6244 out <<
" <Piece Source=\"" << piece_name <<
"\"/>\n";
6246 out <<
" </PUnstructuredGrid>\n";
6247 out <<
"</VTKFile>\n";
6260 const std::vector<std::pair<double, std::string>> ×_and_names)
6264 out <<
"<?xml version=\"1.0\"?>\n";
6267 out <<
"#This file was generated by the deal.II library"
6272 <<
"<VTKFile type=\"Collection\" version=\"0.1\" ByteOrder=\"LittleEndian\">\n";
6273 out <<
" <Collection>\n";
6275 std::streamsize ss = out.precision();
6278 for (
const auto &time_and_name : times_and_names)
6279 out <<
" <DataSet timestep=\"" << time_and_name.first
6280 <<
"\" group=\"\" part=\"0\" file=\"" << time_and_name.second
6283 out <<
" </Collection>\n";
6284 out <<
"</VTKFile>\n";
6296 const std::vector<std::string> &piece_names)
6298 out <<
"!NBLOCKS " << piece_names.size() <<
'\n';
6299 for (
const auto &piece_name : piece_names)
6300 out << piece_name <<
'\n';
6309 const std::vector<std::vector<std::string>> &piece_names)
6313 if (piece_names.empty())
6316 const double nblocks = piece_names[0].size();
6318 ExcMessage(
"piece_names should be a vector of nonempty vectors."));
6320 out <<
"!NBLOCKS " << nblocks <<
'\n';
6321 for (
const auto &domain : piece_names)
6323 Assert(domain.size() == nblocks,
6325 "piece_names should be a vector of equal sized vectors."));
6326 for (
const auto &subdomain : domain)
6327 out << subdomain <<
'\n';
6338 const std::vector<std::pair<
double, std::vector<std::string>>>
6339 ×_and_piece_names)
6343 if (times_and_piece_names.empty())
6346 const double nblocks = times_and_piece_names[0].second.size();