57 #ifdef DEAL_II_WITH_ZLIB 61 #ifdef DEAL_II_WITH_HDF5 75 <<
"Unexpected input: expected line\n <" << arg1
76 <<
">\nbut got\n <" << arg2 <<
">");
82 #ifdef DEAL_II_WITH_ZLIB 88 get_zlib_compression_level(
94 return Z_NO_COMPRESSION;
98 return Z_BEST_COMPRESSION;
100 return Z_DEFAULT_COMPRESSION;
103 return Z_NO_COMPRESSION;
111 template <
typename T>
113 write_compressed_block(
const std::vector<T> & data,
115 std::ostream & output_stream)
117 if (data.size() != 0)
120 auto compressed_data_length = compressBound(data.size() *
sizeof(
T));
121 std::vector<unsigned char> compressed_data(compressed_data_length);
124 compress2(&compressed_data[0],
125 &compressed_data_length,
126 reinterpret_cast<const Bytef *>(data.data()),
127 data.size() *
sizeof(
T),
133 compressed_data.resize(compressed_data_length);
136 const uint32_t compression_header[4] = {
138 static_cast<uint32_t
>(data.size() *
sizeof(
T)),
139 static_cast<uint32_t
>(data.size() *
141 static_cast<uint32_t
>(
142 compressed_data_length)};
144 const auto header_start =
145 reinterpret_cast<const unsigned char *
>(&compression_header[0]);
148 {header_start, header_start + 4 *
sizeof(uint32_t)})
202 return depth > e.depth;
243 return depth > e.depth;
258 template <
int dim,
int spacedim,
typename Number =
double>
260 write_gmv_reorder_data_vectors(
261 const std::vector<Patch<dim, spacedim>> &patches,
265 if (patches.size() == 0)
275 const unsigned int n_data_sets = patches[0].points_are_available ?
276 (patches[0].data.n_rows() - spacedim) :
277 patches[0].data.n_rows();
282 unsigned int next_value = 0;
283 for (
const auto &patch : patches)
285 const unsigned int n_subdivisions = patch.n_subdivisions;
286 (void)n_subdivisions;
288 Assert((patch.data.n_rows() == n_data_sets &&
289 !patch.points_are_available) ||
290 (patch.data.n_rows() == n_data_sets + spacedim &&
291 patch.points_are_available),
293 (n_data_sets + spacedim) :
295 patch.data.n_rows()));
296 Assert(patch.reference_cell != ReferenceCells::get_hypercube<dim>() ||
297 (n_data_sets == 0) ||
298 (patch.data.n_cols() ==
299 Utilities::fixed_power<dim>(n_subdivisions + 1)),
301 n_subdivisions + 1));
303 for (
unsigned int i = 0; i < patch.data.n_cols(); ++i, ++next_value)
304 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
305 data_vectors[data_set][next_value] = patch.data(data_set, i);
308 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
338 for (
unsigned int d = 0;
d < dim; ++
d)
342 unsigned int internal_ind;
353 internal_ind = it->second;
363 const unsigned int pt_index)
382 node_data[node_dim * existing_point.second +
d] =
383 existing_point.first(
d);
391 std::vector<unsigned int> &cell_data)
const 397 cell_data[filtered_cell.first] =
398 filtered_cell.second + local_node_offset;
467 const unsigned int start,
468 const unsigned int d1,
469 const unsigned int d2,
470 const unsigned int d3)
474 const unsigned int base_entry =
500 const unsigned int start,
501 const unsigned int n_points)
505 const unsigned int base_entry = index * n_points;
507 for (
unsigned int i = 0; i < n_points; ++i)
517 const unsigned int dimension,
518 const unsigned int set_num,
521 unsigned int new_dim;
525 if (flags.xdmf_hdf5_output && dimension != 1)
540 for (
unsigned int d = 0;
d < new_dim; ++
d)
543 data_sets.back()[r * new_dim +
d] = data_vectors(set_num +
d, i);
559 const char *gmv_cell_type[4] = {
"",
"line 2",
"quad 4",
"hex 8"};
561 const char *ucd_cell_type[4] = {
"pt",
"line",
"quad",
"hex"};
563 const char *tecplot_cell_type[4] = {
"",
"lineseg",
"quadrilateral",
"brick"};
565 #ifdef DEAL_II_HAVE_TECPLOT 566 const unsigned int tecplot_binary_cell_type[4] = {0, 0, 1, 3};
585 VTK_QUADRATIC_EDGE = 21,
586 VTK_QUADRATIC_TRIANGLE = 22,
587 VTK_QUADRATIC_QUAD = 23,
588 VTK_QUADRATIC_TETRA = 24,
589 VTK_QUADRATIC_HEXAHEDRON = 25,
590 VTK_QUADRATIC_WEDGE = 26,
591 VTK_QUADRATIC_PYRAMID = 27
596 VTK_LAGRANGE_CURVE = 68,
597 VTK_LAGRANGE_TRIANGLE = 69,
598 VTK_LAGRANGE_QUADRILATERAL = 70,
599 VTK_LAGRANGE_TETRAHEDRON = 71,
600 VTK_LAGRANGE_HEXAHEDRON = 72,
601 VTK_LAGRANGE_WEDGE = 73,
602 VTK_LAGRANGE_PYRAMID = 74
609 template <
int dim,
int spacedim>
610 std::array<unsigned int, 3>
612 const bool write_higher_order_cells)
614 std::array<unsigned int, 3> vtk_cell_id{};
616 if (write_higher_order_cells)
620 const std::array<unsigned int, 4> cell_type_by_dim{
623 VTK_LAGRANGE_QUADRILATERAL,
624 VTK_LAGRANGE_HEXAHEDRON}};
625 vtk_cell_id[0] = cell_type_by_dim[dim];
630 vtk_cell_id[0] = VTK_LAGRANGE_TRIANGLE;
639 patch.
data.n_cols() == 3)
641 vtk_cell_id[0] = VTK_TRIANGLE;
645 patch.
data.n_cols() == 6)
647 vtk_cell_id[0] = VTK_QUADRATIC_TRIANGLE;
651 patch.
data.n_cols() == 4)
653 vtk_cell_id[0] = VTK_TETRA;
657 patch.
data.n_cols() == 10)
659 vtk_cell_id[0] = VTK_QUADRATIC_TETRA;
663 patch.
data.n_cols() == 6)
665 vtk_cell_id[0] = VTK_WEDGE;
669 patch.
data.n_cols() == 5)
671 vtk_cell_id[0] = VTK_PYRAMID;
674 else if (patch.
reference_cell == ReferenceCells::get_hypercube<dim>())
676 const std::array<unsigned int, 4> cell_type_by_dim{
677 {VTK_VERTEX, VTK_LINE, VTK_QUAD, VTK_HEXAHEDRON}};
678 vtk_cell_id[0] = cell_type_by_dim[dim];
686 if (patch.
reference_cell != ReferenceCells::get_hypercube<dim>() ||
687 write_higher_order_cells)
688 vtk_cell_id[2] = patch.
data.n_cols();
701 template <
int dim,
int spacedim>
704 const unsigned int xstep,
705 const unsigned int ystep,
706 const unsigned int zstep,
707 const unsigned int n_subdivisions)
712 unsigned int point_no = 0;
717 point_no += (n_subdivisions + 1) * (n_subdivisions + 1) * zstep;
721 point_no += (n_subdivisions + 1) * ystep;
734 for (
unsigned int d = 0;
d < spacedim; ++
d)
735 node[
d] = patch.
data(patch.
data.
size(0) - spacedim +
d, point_no);
744 const double stepsize = 1. / n_subdivisions,
745 xfrac = xstep * stepsize;
751 const double yfrac = ystep * stepsize;
753 node += ((patch.
vertices[3] * xfrac) +
754 (patch.
vertices[2] * (1 - xfrac))) *
758 const double zfrac = zstep * stepsize;
760 node += (((patch.
vertices[5] * xfrac) +
761 (patch.
vertices[4] * (1 - xfrac))) *
764 (patch.
vertices[6] * (1 - xfrac))) *
782 vtk_point_index_from_ijk(
const unsigned i,
785 const std::array<unsigned, 2> &order)
787 const bool ibdy = (i == 0 || i == order[0]);
788 const bool jbdy = (j == 0 || j == order[1]);
790 const int nbdy = (ibdy ? 1 : 0) + (jbdy ? 1 : 0);
794 return (i ? (j ? 2 : 1) : (j ? 3 : 0));
802 return (i - 1) + (j ? order[0] - 1 + order[1] - 1 : 0) + offset;
808 (i ? order[0] - 1 : 2 * (order[0] - 1) + order[1] - 1) +
813 offset += 2 * (order[0] - 1 + order[1] - 1);
815 return offset + (i - 1) + (order[0] - 1) * ((j - 1));
826 vtk_point_index_from_ijk(
const unsigned i,
829 const std::array<unsigned, 3> &order)
831 const bool ibdy = (i == 0 || i == order[0]);
832 const bool jbdy = (j == 0 || j == order[1]);
833 const bool kbdy = (k == 0 || k == order[2]);
835 const int nbdy = (ibdy ? 1 : 0) + (jbdy ? 1 : 0) + (kbdy ? 1 : 0);
839 return (i ? (j ? 2 : 1) : (j ? 3 : 0)) + (k ? 4 : 0);
847 return (i - 1) + (j ? order[0] - 1 + order[1] - 1 : 0) +
848 (k ? 2 * (order[0] - 1 + order[1] - 1) : 0) + offset;
853 (i ? order[0] - 1 : 2 * (order[0] - 1) + order[1] - 1) +
854 (k ? 2 * (order[0] - 1 + order[1] - 1) : 0) + offset;
857 offset += 4 * (order[0] - 1) + 4 * (order[1] - 1);
858 return (k - 1) + (order[2] - 1) * (i ? (j ? 3 : 1) : (j ? 2 : 0)) +
862 offset += 4 * (order[0] - 1 + order[1] - 1 + order[2] - 1);
867 return (j - 1) + ((order[1] - 1) * (k - 1)) +
868 (i ? (order[1] - 1) * (order[2] - 1) : 0) + offset;
870 offset += 2 * (order[1] - 1) * (order[2] - 1);
873 return (i - 1) + ((order[0] - 1) * (k - 1)) +
874 (j ? (order[2] - 1) * (order[0] - 1) : 0) + offset;
876 offset += 2 * (order[2] - 1) * (order[0] - 1);
878 return (i - 1) + ((order[0] - 1) * (j - 1)) +
879 (k ? (order[0] - 1) * (order[1] - 1) : 0) + offset;
884 2 * ((order[1] - 1) * (order[2] - 1) + (order[2] - 1) * (order[0] - 1) +
885 (order[0] - 1) * (order[1] - 1));
886 return offset + (i - 1) +
887 (order[0] - 1) * ((j - 1) + (order[1] - 1) * ((k - 1)));
891 vtk_point_index_from_ijk(
const unsigned,
894 const std::array<unsigned, 0> &)
901 vtk_point_index_from_ijk(
const unsigned,
904 const std::array<unsigned, 1> &)
911 template <
int dim,
int spacedim>
919 for (
const auto &patch : patches)
936 template <
int dim,
int spacedim>
939 const bool write_higher_order_cells,
940 unsigned int &n_nodes,
941 unsigned int &n_cells,
942 unsigned int &n_points_and_n_cells)
946 n_points_and_n_cells = 0;
948 for (
const auto &patch : patches)
955 if (write_higher_order_cells)
958 n_points_and_n_cells +=
964 n_points_and_n_cells +=
971 n_nodes += patch.
data.n_cols();
973 n_points_and_n_cells += patch.
data.n_cols() + 1;
983 template <
typename FlagsType>
990 StreamBase(std::ostream &stream,
const FlagsType &flags)
1005 ExcMessage(
"The derived class you are using needs to " 1006 "reimplement this function if you want to call " 1027 const unsigned int ,
1028 const unsigned int ,
1029 const unsigned int ,
1030 const unsigned int )
1033 ExcMessage(
"The derived class you are using needs to " 1034 "reimplement this function if you want to call " 1046 const unsigned int start,
1047 const unsigned int n_points)
1054 ExcMessage(
"The derived class you are using needs to " 1055 "reimplement this function if you want to call " 1073 template <
typename T>
1087 unsigned int selected_component;
1094 std::ostream &stream;
1099 const FlagsType
flags;
1105 class DXStream :
public StreamBase<DataOutBase::DXFlags>
1125 const unsigned int start,
1126 const unsigned int x_offset,
1127 const unsigned int y_offset,
1128 const unsigned int z_offset);
1136 template <
typename data>
1138 write_dataset(
const unsigned int index,
const std::vector<data> &
values);
1144 class GmvStream :
public StreamBase<DataOutBase::GmvFlags>
1164 const unsigned int start,
1165 const unsigned int x_offset,
1166 const unsigned int y_offset,
1167 const unsigned int z_offset);
1173 class TecplotStream :
public StreamBase<DataOutBase::TecplotFlags>
1193 const unsigned int start,
1194 const unsigned int x_offset,
1195 const unsigned int y_offset,
1196 const unsigned int z_offset);
1202 class UcdStream :
public StreamBase<DataOutBase::UcdFlags>
1224 const unsigned int start,
1225 const unsigned int x_offset,
1226 const unsigned int y_offset,
1227 const unsigned int z_offset);
1235 template <
typename data>
1237 write_dataset(
const unsigned int index,
const std::vector<data> &
values);
1243 class VtkStream :
public StreamBase<DataOutBase::VtkFlags>
1263 const unsigned int start,
1264 const unsigned int x_offset,
1265 const unsigned int y_offset,
1266 const unsigned int z_offset);
1273 const unsigned int start,
1274 const unsigned int n_points);
1285 write_high_order_cell(
const unsigned int index,
1286 const unsigned int start,
1287 const std::vector<unsigned> &connectivity);
1291 class VtuStream :
public StreamBase<DataOutBase::VtkFlags>
1314 const unsigned int start,
1315 const unsigned int x_offset,
1316 const unsigned int y_offset,
1317 const unsigned int z_offset);
1324 const unsigned int start,
1325 const unsigned int n_points);
1336 write_high_order_cell(
const unsigned int index,
1337 const unsigned int start,
1338 const std::vector<unsigned> &connectivity);
1343 template <
typename T>
1354 template <
typename T>
1356 operator<<(const std::vector<T> &);
1368 std::vector<int32_t> cells;
1375 : StreamBase<DataOutBase::DXFlags>(out, f)
1381 DXStream::write_point(
const unsigned int,
const Point<dim> &p)
1383 if (flags.coordinates_binary)
1386 for (
unsigned int d = 0;
d < dim; ++
d)
1388 stream.write(reinterpret_cast<const char *>(data), dim *
sizeof(*data));
1392 for (
unsigned int d = 0;
d < dim; ++
d)
1393 stream << p(
d) <<
'\t';
1402 DXStream::write_cell(
unsigned int,
1408 int nodes[1 << dim];
1412 nodes[GeometryInfo<dim>::dx_to_deal[1]] = start + d1;
1416 nodes[GeometryInfo<dim>::dx_to_deal[2]] = start + d2;
1417 nodes[GeometryInfo<dim>::dx_to_deal[3]] = start + d2 + d1;
1421 nodes[GeometryInfo<dim>::dx_to_deal[4]] = start + d3;
1422 nodes[GeometryInfo<dim>::dx_to_deal[5]] = start + d3 + d1;
1423 nodes[GeometryInfo<dim>::dx_to_deal[6]] = start + d3 + d2;
1424 nodes[GeometryInfo<dim>::dx_to_deal[7]] = start + d3 + d2 + d1;
1429 if (flags.int_binary)
1430 stream.write(reinterpret_cast<const char *>(nodes),
1431 (1 << dim) *
sizeof(*nodes));
1434 const unsigned int final = (1 << dim) - 1;
1435 for (
unsigned int i = 0; i <
final; ++i)
1436 stream << nodes[i] <<
'\t';
1437 stream << nodes[
final] <<
'\n';
1443 template <
typename data>
1445 DXStream::write_dataset(
const unsigned int,
const std::vector<data> &
values)
1447 if (flags.data_binary)
1449 stream.write(reinterpret_cast<const char *>(values.data()),
1450 values.size() *
sizeof(data));
1454 for (
unsigned int i = 0; i < values.size(); ++i)
1455 stream <<
'\t' << values[i];
1465 : StreamBase<DataOutBase::GmvFlags>(out, f)
1471 GmvStream::write_point(
const unsigned int,
const Point<dim> &p)
1475 stream << p(selected_component) <<
' ';
1482 GmvStream::write_cell(
unsigned int,
1489 const unsigned int start = s + 1;
1490 stream << gmv_cell_type[dim] <<
'\n';
1495 stream <<
'\t' << start + d1;
1498 stream <<
'\t' << start + d2 + d1 <<
'\t' << start + d2;
1501 stream <<
'\t' << start + d3 <<
'\t' << start + d3 + d1 <<
'\t' 1502 << start + d3 + d2 + d1 <<
'\t' << start + d3 + d2;
1511 TecplotStream::TecplotStream(std::ostream & out,
1513 : StreamBase<DataOutBase::TecplotFlags>(out, f)
1519 TecplotStream::write_point(
const unsigned int,
const Point<dim> &p)
1523 stream << p(selected_component) <<
'\n';
1530 TecplotStream::write_cell(
unsigned int,
1536 const unsigned int start = s + 1;
1541 stream <<
'\t' << start + d1;
1544 stream <<
'\t' << start + d2 + d1 <<
'\t' << start + d2;
1547 stream <<
'\t' << start + d3 <<
'\t' << start + d3 + d1 <<
'\t' 1548 << start + d3 + d2 + d1 <<
'\t' << start + d3 + d2;
1558 : StreamBase<DataOutBase::UcdFlags>(out, f)
1564 UcdStream::write_point(
const unsigned int index,
const Point<dim> &p)
1566 stream << index + 1 <<
" ";
1568 for (
unsigned int i = 0; i < dim; ++i)
1569 stream << p(i) <<
' ';
1571 for (
unsigned int i = dim; i < 3; ++i)
1580 UcdStream::write_cell(
unsigned int index,
1586 int nodes[1 << dim];
1590 nodes[GeometryInfo<dim>::ucd_to_deal[1]] = start + d1;
1594 nodes[GeometryInfo<dim>::ucd_to_deal[2]] = start + d2;
1595 nodes[GeometryInfo<dim>::ucd_to_deal[3]] = start + d2 + d1;
1599 nodes[GeometryInfo<dim>::ucd_to_deal[4]] = start + d3;
1600 nodes[GeometryInfo<dim>::ucd_to_deal[5]] = start + d3 + d1;
1601 nodes[GeometryInfo<dim>::ucd_to_deal[6]] = start + d3 + d2;
1602 nodes[GeometryInfo<dim>::ucd_to_deal[7]] = start + d3 + d2 + d1;
1608 stream << index + 1 <<
"\t0 " << ucd_cell_type[dim];
1609 const unsigned int final = (1 << dim);
1610 for (
unsigned int i = 0; i <
final; ++i)
1611 stream <<
'\t' << nodes[i] + 1;
1617 template <
typename data>
1619 UcdStream::write_dataset(
const unsigned int index,
1620 const std::vector<data> &values)
1622 stream << index + 1;
1623 for (
unsigned int i = 0; i < values.size(); ++i)
1624 stream <<
'\t' << values[i];
1633 : StreamBase<DataOutBase::VtkFlags>(out, f)
1639 VtkStream::write_point(
const unsigned int,
const Point<dim> &p)
1644 for (
unsigned int i = dim; i < 3; ++i)
1653 VtkStream::write_cell(
unsigned int,
1659 stream << GeometryInfo<dim>::vertices_per_cell <<
'\t' << start;
1662 stream <<
'\t' << start + d1;
1666 stream <<
'\t' << start + d2 + d1 <<
'\t' << start + d2;
1669 stream <<
'\t' << start + d3 <<
'\t' << start + d3 + d1 <<
'\t' 1670 << start + d3 + d2 + d1 <<
'\t' << start + d3 + d2;
1678 VtkStream::write_cell_single(
const unsigned int index,
1679 const unsigned int start,
1680 const unsigned int n_points)
1684 stream <<
'\t' << n_points;
1685 for (
unsigned int i = 0; i < n_points; ++i)
1686 stream <<
'\t' << start + i;
1692 VtkStream::write_high_order_cell(
const unsigned int,
1693 const unsigned int start,
1694 const std::vector<unsigned> &connectivity)
1696 stream << connectivity.size();
1697 for (
const auto &c : connectivity)
1698 stream <<
'\t' << start + c;
1703 : StreamBase<DataOutBase::VtkFlags>(out, f)
1709 VtuStream::write_point(
const unsigned int,
const Point<dim> &p)
1711 #if !defined(DEAL_II_WITH_ZLIB) 1715 for (
unsigned int i = dim; i < 3; ++i)
1720 for (
unsigned int i = 0; i < dim; ++i)
1722 for (
unsigned int i = dim; i < 3; ++i)
1729 VtuStream::flush_points()
1731 #ifdef DEAL_II_WITH_ZLIB 1742 VtuStream::write_cell(
unsigned int,
1748 #if !defined(DEAL_II_WITH_ZLIB) 1752 stream <<
'\t' << start + d1;
1755 stream <<
'\t' << start + d2 + d1 <<
'\t' << start + d2;
1758 stream <<
'\t' << start + d3 <<
'\t' << start + d3 + d1 <<
'\t' 1759 << start + d3 + d2 + d1 <<
'\t' << start + d3 + d2;
1765 cells.push_back(start);
1768 cells.push_back(start + d1);
1771 cells.push_back(start + d2 + d1);
1772 cells.push_back(start + d2);
1775 cells.push_back(start + d3);
1776 cells.push_back(start + d3 + d1);
1777 cells.push_back(start + d3 + d2 + d1);
1778 cells.push_back(start + d3 + d2);
1786 VtuStream::write_cell_single(
const unsigned int index,
1787 const unsigned int start,
1788 const unsigned int n_points)
1792 #if !defined(DEAL_II_WITH_ZLIB) 1793 for (
unsigned int i = 0; i < n_points; ++i)
1794 stream <<
'\t' << start + i;
1797 for (
unsigned int i = 0; i < n_points; ++i)
1798 cells.push_back(start + i);
1804 VtuStream::write_high_order_cell(
const unsigned int,
1805 const unsigned int start,
1806 const std::vector<unsigned> &connectivity)
1808 #if !defined(DEAL_II_WITH_ZLIB) 1809 for (
const auto &c : connectivity)
1810 stream <<
'\t' << start + c;
1813 for (
const auto &c : connectivity)
1814 cells.push_back(start + c);
1819 VtuStream::flush_cells()
1821 #ifdef DEAL_II_WITH_ZLIB 1824 *
this << cells <<
'\n';
1830 template <
typename T>
1832 VtuStream::operator<<(const std::vector<T> &data)
1834 #ifdef DEAL_II_WITH_ZLIB 1837 write_compressed_block(data, flags, stream);
1839 for (
unsigned int i = 0; i < data.size(); ++i)
1840 stream << data[i] <<
' ';
1854 template <
int dim,
int spacedim>
1858 template <
int dim,
int spacedim>
1862 template <
int dim,
int spacedim>
1864 : patch_index(no_neighbor)
1866 , points_are_available(false)
1872 neighbors[i] = no_neighbor;
1880 template <
int dim,
int spacedim>
1903 if (data.n_rows() != patch.
data.n_rows())
1906 if (data.n_cols() != patch.
data.n_cols())
1909 for (
unsigned int i = 0; i < data.n_rows(); ++i)
1910 for (
unsigned int j = 0; j < data.n_cols(); ++j)
1911 if (data[i][j] != patch.
data[i][j])
1919 template <
int dim,
int spacedim>
1925 sizeof(neighbors) /
sizeof(neighbors[0]) *
1935 template <
int dim,
int spacedim>
1943 data.swap(other_patch.
data);
1950 template <
int spacedim>
1954 template <
int spacedim>
1958 template <
int spacedim>
1960 Patch<0, spacedim>::no_neighbor};
1962 template <
int spacedim>
1965 template <
int spacedim>
1967 : patch_index(no_neighbor)
1968 , points_are_available(false)
1976 template <
int spacedim>
1980 const unsigned int dim = 0;
1994 if (data.n_rows() != patch.
data.n_rows())
1997 if (data.n_cols() != patch.
data.n_cols())
2000 for (
unsigned int i = 0; i < data.n_rows(); ++i)
2001 for (
unsigned int j = 0; j < data.n_cols(); ++j)
2002 if (data[i][j] != patch.
data[i][j])
2010 template <
int spacedim>
2022 template <
int spacedim>
2027 data.swap(other_patch.
data);
2034 : write_preamble(write_preamble)
2041 space_dimension_labels.emplace_back(
"x");
2042 space_dimension_labels.emplace_back(
"y");
2043 space_dimension_labels.emplace_back(
"z");
2049 : space_dimension_labels(labels)
2063 const bool bicubic_patch,
2064 const bool external_data)
2066 , bicubic_patch(bicubic_patch)
2067 , external_data(external_data)
2072 const bool xdmf_hdf5_output)
2073 : filter_duplicate_vertices(filter_duplicate_vertices)
2074 , xdmf_hdf5_output(xdmf_hdf5_output)
2082 "Filter duplicate vertices",
2085 "Whether to remove duplicate vertex values. deal.II duplicates " 2086 "vertices once for each adjacent cell so that it can output " 2087 "discontinuous quantities for which there may be more than one " 2088 "value for each vertex position. Setting this flag to " 2089 "'true' will merge all of these values by selecting a " 2090 "random one and outputting this as 'the' value for the vertex. " 2091 "As long as the data to be output corresponds to continuous " 2092 "fields, merging vertices has no effect. On the other hand, " 2093 "if the data to be output corresponds to discontinuous fields " 2094 "(either because you are using a discontinuous finite element, " 2095 "or because you are using a DataPostprocessor that yields " 2096 "discontinuous data, or because the data to be output has been " 2097 "produced by entirely different means), then the data in the " 2098 "output file no longer faithfully represents the underlying data " 2099 "because the discontinuous field has been replaced by a " 2100 "continuous one. Note also that the filtering can not occur " 2101 "on processor boundaries. Thus, a filtered discontinuous field " 2102 "looks like a continuous field inside of a subdomain, " 2103 "but like a discontinuous field at the subdomain boundary." 2105 "In any case, filtering results in drastically smaller output " 2106 "files (smaller by about a factor of 2^dim).");
2111 "Whether the data will be used in an XDMF/HDF5 combination.");
2126 const bool int_binary,
2127 const bool coordinates_binary,
2128 const bool data_binary)
2129 : write_neighbors(write_neighbors)
2130 , int_binary(int_binary)
2131 , coordinates_binary(coordinates_binary)
2132 , data_binary(data_binary)
2133 , data_double(false)
2143 "A boolean field indicating whether neighborship " 2144 "information between cells is to be written to the " 2145 "OpenDX output file");
2149 "Output format of integer numbers, which is " 2150 "either a text representation (ascii) or binary integer " 2151 "values of 32 or 64 bits length");
2155 "Output format of vertex coordinates, which is " 2156 "either a text representation (ascii) or binary " 2157 "floating point values of 32 or 64 bits length");
2161 "Output format of data values, which is " 2162 "either a text representation (ascii) or binary " 2163 "floating point values of 32 or 64 bits length");
2183 "A flag indicating whether a comment should be " 2184 "written to the beginning of the output file " 2185 "indicating date and time of creation as well " 2186 "as the creating program");
2194 write_preamble = prm.
get_bool(
"Write preamble");
2200 const int azimuth_angle,
2201 const int polar_angle,
2202 const unsigned int line_thickness,
2204 const bool draw_colorbar)
2207 , height_vector(height_vector)
2208 , azimuth_angle(azimuth_angle)
2209 , polar_angle(polar_angle)
2210 , line_thickness(line_thickness)
2212 , draw_colorbar(draw_colorbar)
2223 "A flag indicating whether POVRAY should use smoothed " 2224 "triangles instead of the usual ones");
2228 "Whether POVRAY should use bicubic patches");
2232 "Whether camera and lighting information should " 2233 "be put into an external file \"data.inc\" or into " 2234 "the POVRAY input file");
2242 smooth = prm.
get_bool(
"Use smooth triangles");
2243 bicubic_patch = prm.
get_bool(
"Use bicubic patches");
2244 external_data = prm.
get_bool(
"Include external file");
2250 const unsigned int color_vector,
2252 const unsigned int size,
2253 const double line_width,
2254 const double azimut_angle,
2255 const double turn_angle,
2256 const double z_scaling,
2257 const bool draw_mesh,
2258 const bool draw_cells,
2259 const bool shade_cells,
2261 : height_vector(height_vector)
2262 , color_vector(color_vector)
2263 , size_type(size_type)
2265 , line_width(line_width)
2266 , azimut_angle(azimut_angle)
2267 , turn_angle(turn_angle)
2268 , z_scaling(z_scaling)
2269 , draw_mesh(draw_mesh)
2270 , draw_cells(draw_cells)
2271 , shade_cells(shade_cells)
2272 , color_function(color_function)
2311 double sum = xmax + xmin;
2312 double sum13 = xmin + 3 * xmax;
2313 double sum22 = 2 * xmin + 2 * xmax;
2314 double sum31 = 3 * xmin + xmax;
2315 double dif = xmax - xmin;
2316 double rezdif = 1.0 / dif;
2320 if (x < (sum31) / 4)
2322 else if (x < (sum22) / 4)
2324 else if (x < (sum13) / 4)
2335 rgb_values.
green = 0;
2336 rgb_values.
blue = (x - xmin) * 4. * rezdif;
2340 rgb_values.
green = (4 * x - 3 * xmin - xmax) * rezdif;
2341 rgb_values.
blue = (sum22 - 4. * x) * rezdif;
2344 rgb_values.
red = (4 * x - 2 *
sum) * rezdif;
2345 rgb_values.
green = (xmin + 3 * xmax - 4 * x) * rezdif;
2346 rgb_values.
blue = 0;
2350 rgb_values.
green = (4 * x - xmin - 3 * xmax) * rezdif;
2351 rgb_values.
blue = (4. * x - sum13) * rezdif;
2358 rgb_values.
red = rgb_values.
green = rgb_values.
blue = 1;
2372 (x - xmin) / (xmax - xmin);
2385 1 - (x - xmin) / (xmax - xmin);
2397 "Number of the input vector that is to be used to " 2398 "generate height information");
2402 "Number of the input vector that is to be used to " 2403 "generate color information");
2407 "Whether width or height should be scaled to match " 2412 "The size (width or height) to which the eps output " 2413 "file is to be scaled");
2417 "The width in which the postscript renderer is to " 2422 "Angle of the viewing position against the vertical " 2427 "Angle of the viewing direction against the y-axis");
2431 "Scaling for the z-direction relative to the scaling " 2432 "used in x- and y-directions");
2436 "Whether the mesh lines, or only the surface should be " 2441 "Whether only the mesh lines, or also the interior of " 2442 "cells should be plotted. If this flag is false, then " 2443 "one can see through the mesh");
2447 "Whether the interior of cells shall be shaded");
2451 "default|grey scale|reverse grey scale"),
2452 "Name of a color function used to colorize mesh lines " 2453 "and/or cell interiors");
2463 if (prm.
get(
"Scale to width or height") ==
"width")
2475 if (prm.
get(
"Color function") ==
"default")
2477 else if (prm.
get(
"Color function") ==
"grey scale")
2479 else if (prm.
get(
"Color function") ==
"reverse grey scale")
2490 : zone_name(zone_name)
2491 , solution_time(solution_time)
2505 const unsigned int cycle,
2506 const bool print_date_and_time,
2508 const bool write_higher_order_cells)
2511 , print_date_and_time(print_date_and_time)
2512 , compression_level(compression_level)
2513 , write_higher_order_cells(write_higher_order_cells)
2521 if (format_name ==
"none")
2524 if (format_name ==
"dx")
2527 if (format_name ==
"ucd")
2530 if (format_name ==
"gnuplot")
2533 if (format_name ==
"povray")
2536 if (format_name ==
"eps")
2539 if (format_name ==
"gmv")
2542 if (format_name ==
"tecplot")
2545 if (format_name ==
"tecplot_binary")
2548 if (format_name ==
"vtk")
2551 if (format_name ==
"vtu")
2554 if (format_name ==
"deal.II intermediate")
2557 if (format_name ==
"hdf5")
2561 ExcMessage(
"The given file format name is not recognized: <" +
2562 format_name +
">"));
2573 return "none|dx|ucd|gnuplot|povray|eps|gmv|tecplot|tecplot_binary|vtk|vtu|hdf5|svg|deal.II intermediate";
2581 switch (output_format)
2620 template <
int dim,
int spacedim,
typename StreamType>
2625 unsigned int count = 0;
2627 for (
const auto &patch : patches)
2632 if (patch.reference_cell != ReferenceCells::get_hypercube<dim>())
2636 for (
unsigned int point_no = 0; point_no < patch.data.n_cols();
2639 for (
unsigned int d = 0;
d < spacedim; ++
d)
2641 patch.data(patch.data.size(0) - spacedim +
d, point_no);
2643 out.write_point(count++, node);
2648 const unsigned int n_subdivisions = patch.n_subdivisions;
2649 const unsigned int n = n_subdivisions + 1;
2652 const unsigned int n1 = (dim > 0) ? n : 1;
2653 const unsigned int n2 = (dim > 1) ? n : 1;
2654 const unsigned int n3 = (dim > 2) ? n : 1;
2656 for (
unsigned int i3 = 0; i3 < n3; ++i3)
2657 for (
unsigned int i2 = 0; i2 < n2; ++i2)
2658 for (
unsigned int i1 = 0; i1 < n1; ++i1)
2660 count++, compute_node(patch, i1, i2, i3, n_subdivisions));
2666 template <
int dim,
int spacedim,
typename StreamType>
2671 unsigned int count = 0;
2672 unsigned int first_vertex_of_patch = 0;
2673 for (
const auto &patch : patches)
2676 if (patch.reference_cell != ReferenceCells::get_hypercube<dim>())
2678 out.write_cell_single(count++,
2679 first_vertex_of_patch,
2680 patch.data.n_cols());
2681 first_vertex_of_patch += patch.data.n_cols();
2685 const unsigned int n_subdivisions = patch.n_subdivisions;
2686 const unsigned int n = n_subdivisions + 1;
2688 const unsigned int n1 = (dim > 0) ? n_subdivisions : 1;
2689 const unsigned int n2 = (dim > 1) ? n_subdivisions : 1;
2690 const unsigned int n3 = (dim > 2) ? n_subdivisions : 1;
2692 const unsigned int d1 = 1;
2693 const unsigned int d2 = n;
2694 const unsigned int d3 = n * n;
2695 for (
unsigned int i3 = 0; i3 < n3; ++i3)
2696 for (
unsigned int i2 = 0; i2 < n2; ++i2)
2697 for (
unsigned int i1 = 0; i1 < n1; ++i1)
2699 const unsigned int offset =
2700 first_vertex_of_patch + i3 * d3 + i2 * d2 + i1 * d1;
2702 out.template write_cell<dim>(count++, offset, d1, d2, d3);
2705 first_vertex_of_patch +=
2706 Utilities::fixed_power<dim>(n_subdivisions + 1);
2713 template <
int dim,
int spacedim,
typename StreamType>
2719 unsigned int first_vertex_of_patch = 0;
2720 unsigned int count = 0;
2722 std::vector<unsigned> connectivity;
2724 std::array<unsigned, dim> cell_order;
2726 for (
const auto &patch : patches)
2728 if (patch.reference_cell != ReferenceCells::get_hypercube<dim>())
2730 connectivity.resize(patch.data.n_cols());
2732 for (
unsigned int i = 0; i < patch.data.n_cols(); ++i)
2733 connectivity[i] = i;
2735 out.template write_high_order_cell<dim>(count++,
2736 first_vertex_of_patch,
2739 first_vertex_of_patch += patch.data.n_cols();
2743 const unsigned int n_subdivisions = patch.n_subdivisions;
2744 const unsigned int n = n_subdivisions + 1;
2746 cell_order.fill(n_subdivisions);
2747 connectivity.resize(Utilities::fixed_power<dim>(n));
2750 const unsigned int n1 = (dim > 0) ? n_subdivisions : 0;
2751 const unsigned int n2 = (dim > 1) ? n_subdivisions : 0;
2752 const unsigned int n3 = (dim > 2) ? n_subdivisions : 0;
2754 const unsigned int d1 = 1;
2755 const unsigned int d2 = n;
2756 const unsigned int d3 = n * n;
2757 for (
unsigned int i3 = 0; i3 <= n3; ++i3)
2758 for (
unsigned int i2 = 0; i2 <= n2; ++i2)
2759 for (
unsigned int i1 = 0; i1 <= n1; ++i1)
2761 const unsigned int local_index =
2762 i3 * d3 + i2 * d2 + i1 * d1;
2763 const unsigned int connectivity_index =
2764 vtk_point_index_from_ijk(i1, i2, i3, cell_order);
2765 connectivity[connectivity_index] = local_index;
2768 out.template write_high_order_cell<dim>(count++,
2769 first_vertex_of_patch,
2773 first_vertex_of_patch += Utilities::fixed_power<dim>(n);
2781 template <
int dim,
int spacedim,
class StreamType>
2784 unsigned int n_data_sets,
2785 const bool double_precision,
2789 unsigned int count = 0;
2791 for (
const auto &patch : patches)
2793 const unsigned int n_subdivisions = patch.n_subdivisions;
2794 const unsigned int n = n_subdivisions + 1;
2796 Assert((patch.data.n_rows() == n_data_sets &&
2797 !patch.points_are_available) ||
2798 (patch.data.n_rows() == n_data_sets + spacedim &&
2799 patch.points_are_available),
2801 (n_data_sets + spacedim) :
2803 patch.data.n_rows()));
2804 Assert(patch.data.n_cols() == Utilities::fixed_power<dim>(n),
2807 std::vector<float> floats(n_data_sets);
2808 std::vector<double> doubles(n_data_sets);
2811 for (
unsigned int i = 0; i < Utilities::fixed_power<dim>(n);
2813 if (double_precision)
2815 for (
unsigned int data_set = 0; data_set < n_data_sets;
2817 doubles[data_set] = patch.data(data_set, i);
2818 out.write_dataset(count, doubles);
2822 for (
unsigned int data_set = 0; data_set < n_data_sets;
2824 floats[data_set] = patch.data(data_set, i);
2825 out.write_dataset(count, floats);
2849 camera_vertical[0] = camera_horizontal[1] * camera_direction[2] -
2850 camera_horizontal[2] * camera_direction[1];
2851 camera_vertical[1] = camera_horizontal[2] * camera_direction[0] -
2852 camera_horizontal[0] * camera_direction[2];
2853 camera_vertical[2] = camera_horizontal[0] * camera_direction[1] -
2854 camera_horizontal[1] * camera_direction[0];
2858 phi /= (point[0] - camera_position[0]) * camera_direction[0] +
2859 (point[1] - camera_position[1]) * camera_direction[1] +
2860 (point[2] - camera_position[2]) * camera_direction[2];
2864 camera_position[0] + phi * (point[0] - camera_position[0]);
2866 camera_position[1] + phi * (point[1] - camera_position[1]);
2868 camera_position[2] + phi * (point[2] - camera_position[2]);
2871 projection_decomposition[0] = (projection[0] - camera_position[0] -
2872 camera_focus * camera_direction[0]) *
2873 camera_horizontal[0];
2874 projection_decomposition[0] += (projection[1] - camera_position[1] -
2875 camera_focus * camera_direction[1]) *
2876 camera_horizontal[1];
2877 projection_decomposition[0] += (projection[2] - camera_position[2] -
2878 camera_focus * camera_direction[2]) *
2879 camera_horizontal[2];
2881 projection_decomposition[1] = (projection[0] - camera_position[0] -
2882 camera_focus * camera_direction[0]) *
2884 projection_decomposition[1] += (projection[1] - camera_position[1] -
2885 camera_focus * camera_direction[1]) *
2887 projection_decomposition[1] += (projection[2] - camera_position[2] -
2888 camera_focus * camera_direction[2]) *
2891 return projection_decomposition;
2905 for (
int i = 0; i < 2; ++i)
2907 for (
int j = 0; j < 2 - i; ++j)
2909 if (points[j][2] > points[j + 1][2])
2912 points[j] = points[j + 1];
2913 points[j + 1] = temp;
2920 v_inter = points[1];
2927 A[0][0] = v_max[0] - v_min[0];
2928 A[0][1] = v_inter[0] - v_min[0];
2929 A[1][0] = v_max[1] - v_min[1];
2930 A[1][1] = v_inter[1] - v_min[1];
2936 bool col_change =
false;
2945 double temp = A[1][0];
2950 for (
unsigned int k = 0; k < 1; k++)
2952 for (
unsigned int i = k + 1; i < 2; i++)
2954 x = A[i][k] / A[k][k];
2956 for (
unsigned int j = k + 1; j < 2; j++)
2957 A[i][j] = A[i][j] - A[k][j] * x;
2959 b[i] = b[i] - b[k] * x;
2963 b[1] = b[1] / A[1][1];
2965 for (
int i = 0; i >= 0; i--)
2969 for (
unsigned int j = i + 1; j < 2; j++)
2970 sum = sum - A[i][j] * b[j];
2972 b[i] = sum / A[i][i];
2982 double c = b[0] * (v_max[2] - v_min[2]) + b[1] * (v_inter[2] - v_min[2]) +
2986 A[0][0] = v_max[0] - v_min[0];
2987 A[0][1] = v_inter[0] - v_min[0];
2988 A[1][0] = v_max[1] - v_min[1];
2989 A[1][1] = v_inter[1] - v_min[1];
2991 b[0] = 1.0 - v_min[0];
3003 double temp = A[1][0];
3008 for (
unsigned int k = 0; k < 1; k++)
3010 for (
unsigned int i = k + 1; i < 2; i++)
3012 x = A[i][k] / A[k][k];
3014 for (
unsigned int j = k + 1; j < 2; j++)
3015 A[i][j] = A[i][j] - A[k][j] * x;
3017 b[i] = b[i] - b[k] * x;
3021 b[1] = b[1] / A[1][1];
3023 for (
int i = 0; i >= 0; i--)
3027 for (
unsigned int j = i + 1; j < 2; j++)
3028 sum = sum - A[i][j] * b[j];
3030 b[i] = sum / A[i][i];
3040 gradient[0] = b[0] * (v_max[2] - v_min[2]) +
3041 b[1] * (v_inter[2] - v_min[2]) - c + v_min[2];
3044 A[0][0] = v_max[0] - v_min[0];
3045 A[0][1] = v_inter[0] - v_min[0];
3046 A[1][0] = v_max[1] - v_min[1];
3047 A[1][1] = v_inter[1] - v_min[1];
3050 b[1] = 1.0 - v_min[1];
3061 double temp = A[1][0];
3066 for (
unsigned int k = 0; k < 1; k++)
3068 for (
unsigned int i = k + 1; i < 2; i++)
3070 x = A[i][k] / A[k][k];
3072 for (
unsigned int j = k + 1; j < 2; j++)
3073 A[i][j] = A[i][j] - A[k][j] * x;
3075 b[i] = b[i] - b[k] * x;
3079 b[1] = b[1] / A[1][1];
3081 for (
int i = 0; i >= 0; i--)
3085 for (
unsigned int j = i + 1; j < 2; j++)
3086 sum = sum - A[i][j] * b[j];
3088 b[i] = sum / A[i][i];
3098 gradient[1] = b[0] * (v_max[2] - v_min[2]) +
3099 b[1] * (v_inter[2] - v_min[2]) - c + v_min[2];
3102 double gradient_norm =
3103 std::sqrt(std::pow(gradient[0], 2.0) + std::pow(gradient[1], 2.0));
3104 gradient[0] /= gradient_norm;
3105 gradient[1] /= gradient_norm;
3107 double lambda = -gradient[0] * (v_min[0] - v_max[0]) -
3108 gradient[1] * (v_min[1] - v_max[1]);
3112 gradient_parameters[0] = v_min[0];
3113 gradient_parameters[1] = v_min[1];
3115 gradient_parameters[2] = v_min[0] + lambda * gradient[0];
3116 gradient_parameters[3] = v_min[1] + lambda * gradient[1];
3118 gradient_parameters[4] = v_min[2];
3119 gradient_parameters[5] = v_max[2];
3121 return gradient_parameters;
3127 template <
int dim,
int spacedim>
3131 const std::vector<std::string> & data_names,
3133 std::tuple<
unsigned int,
3146 #ifndef DEAL_II_WITH_MPI 3155 if (patches.size() == 0)
3159 const unsigned int n_data_sets = data_names.size();
3161 UcdStream ucd_out(out, flags);
3164 unsigned int n_nodes;
3166 compute_sizes<dim, spacedim>(patches, n_nodes,
n_cells);
3172 <<
"# This file was generated by the deal.II library." <<
'\n' 3176 <<
"# For a description of the UCD format see the AVS Developer's guide." 3182 out << n_nodes <<
' ' << n_cells <<
' ' << n_data_sets <<
' ' << 0
3195 if (n_data_sets != 0)
3197 out << n_data_sets <<
" ";
3198 for (
unsigned int i = 0; i < n_data_sets; ++i)
3203 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
3204 out << data_names[data_set]
3208 write_data(patches, n_data_sets,
true, ucd_out);
3218 template <
int dim,
int spacedim>
3222 const std::vector<std::string> & data_names,
3224 std::tuple<
unsigned int,
3236 #ifndef DEAL_II_WITH_MPI 3245 if (patches.size() == 0)
3249 DXStream dx_out(out, flags);
3252 unsigned int offset = 0;
3254 const unsigned int n_data_sets = data_names.size();
3257 unsigned int n_nodes;
3259 compute_sizes<dim, spacedim>(patches, n_nodes,
n_cells);
3261 out <<
"object \"vertices\" class array type float rank 1 shape " 3262 << spacedim <<
" items " << n_nodes;
3266 out <<
" lsb ieee data 0" <<
'\n';
3267 offset += n_nodes * spacedim *
sizeof(float);
3271 out <<
" data follows" <<
'\n';
3280 out <<
"object \"cells\" class array type int rank 1 shape " 3285 out <<
" lsb binary data " << offset <<
'\n';
3286 offset += n_cells *
sizeof(int);
3290 out <<
" data follows" <<
'\n';
3296 out <<
"attribute \"element type\" string \"";
3303 out <<
"\"" <<
'\n' <<
"attribute \"ref\" string \"positions\"" <<
'\n';
3310 out <<
"object \"neighbors\" class array type int rank 1 shape " 3314 for (
const auto &patch : patches)
3316 const unsigned int n = patch.n_subdivisions;
3317 const unsigned int n1 = (dim > 0) ? n : 1;
3318 const unsigned int n2 = (dim > 1) ? n : 1;
3319 const unsigned int n3 = (dim > 2) ? n : 1;
3320 unsigned int cells_per_patch = Utilities::fixed_power<dim>(n);
3321 unsigned int dx = 1;
3322 unsigned int dy = n;
3323 unsigned int dz = n * n;
3325 const unsigned int patch_start =
3326 patch.patch_index * cells_per_patch;
3328 for (
unsigned int i3 = 0; i3 < n3; ++i3)
3329 for (
unsigned int i2 = 0; i2 < n2; ++i2)
3330 for (
unsigned int i1 = 0; i1 < n1; ++i1)
3332 const unsigned int nx = i1 *
dx;
3333 const unsigned int ny = i2 * dy;
3334 const unsigned int nz = i3 * dz;
3346 const unsigned int nn = patch.neighbors[0];
3348 if (nn != patch.no_neighbor)
3350 << (nn * cells_per_patch + ny + nz + dx * (n - 1));
3356 out <<
'\t' << patch_start + nx - dx + ny + nz;
3361 const unsigned int nn = patch.neighbors[1];
3363 if (nn != patch.no_neighbor)
3364 out << (nn * cells_per_patch + ny + nz);
3370 out <<
'\t' << patch_start + nx + dx + ny + nz;
3377 const unsigned int nn = patch.neighbors[2];
3379 if (nn != patch.no_neighbor)
3381 << (nn * cells_per_patch + nx + nz + dy * (n - 1));
3387 out <<
'\t' << patch_start + nx + ny - dy + nz;
3392 const unsigned int nn = patch.neighbors[3];
3394 if (nn != patch.no_neighbor)
3395 out << (nn * cells_per_patch + nx + nz);
3401 out <<
'\t' << patch_start + nx + ny + dy + nz;
3409 const unsigned int nn = patch.neighbors[4];
3411 if (nn != patch.no_neighbor)
3413 << (nn * cells_per_patch + nx + ny + dz * (n - 1));
3419 out <<
'\t' << patch_start + nx + ny + nz - dz;
3424 const unsigned int nn = patch.neighbors[5];
3426 if (nn != patch.no_neighbor)
3427 out << (nn * cells_per_patch + nx + ny);
3433 out <<
'\t' << patch_start + nx + ny + nz + dz;
3441 if (n_data_sets != 0)
3443 out <<
"object \"data\" class array type float rank 1 shape " 3444 << n_data_sets <<
" items " << n_nodes;
3448 out <<
" lsb ieee data " << offset <<
'\n';
3449 offset += n_data_sets * n_nodes *
3450 ((flags.
data_double) ?
sizeof(
double) :
sizeof(float));
3454 out <<
" data follows" <<
'\n';
3459 out <<
"attribute \"dep\" string \"positions\"" <<
'\n';
3463 out <<
"object \"data\" class constantarray type float rank 0 items " 3464 << n_nodes <<
" data follows" <<
'\n' 3470 out <<
"object \"deal data\" class field" <<
'\n' 3471 <<
"component \"positions\" value \"vertices\"" <<
'\n' 3472 <<
"component \"connections\" value \"cells\"" <<
'\n' 3473 <<
"component \"data\" value \"data\"" <<
'\n';
3476 out <<
"component \"neighbors\" value \"neighbors\"" <<
'\n';
3483 out <<
"end" <<
'\n';
3501 template <
int dim,
int spacedim>
3505 const std::vector<std::string> & data_names,
3507 std::tuple<
unsigned int,
3516 #ifndef DEAL_II_WITH_MPI 3526 if (patches.size() == 0)
3530 const unsigned int n_data_sets = data_names.size();
3534 out <<
"# This file was generated by the deal.II library." <<
'\n' 3538 <<
"# For a description of the GNUPLOT format see the GNUPLOT manual." 3545 for (
unsigned int spacedim_n = 0; spacedim_n < spacedim; ++spacedim_n)
3550 for (
const auto &data_name : data_names)
3551 out <<
'<' << data_name <<
"> ";
3557 for (
const auto &patch : patches)
3559 const unsigned int n_subdivisions = patch.n_subdivisions;
3560 const unsigned int n = n_subdivisions + 1;
3562 const unsigned int n1 = (dim > 0) ? n : 1;
3563 const unsigned int n2 = (dim > 1) ? n : 1;
3564 const unsigned int n3 = (dim > 2) ? n : 1;
3565 unsigned int d1 = 1;
3566 unsigned int d2 = n;
3567 unsigned int d3 = n * n;
3569 Assert((patch.data.n_rows() == n_data_sets &&
3570 !patch.points_are_available) ||
3571 (patch.data.n_rows() == n_data_sets + spacedim &&
3572 patch.points_are_available),
3574 (n_data_sets + spacedim) :
3576 patch.data.n_rows()));
3577 Assert(patch.data.n_cols() == Utilities::fixed_power<dim>(n),
3583 for (
unsigned int i2 = 0; i2 < n2; ++i2)
3585 for (
unsigned int i1 = 0; i1 < n1; ++i1)
3588 out << compute_node(patch, i1, i2, 0, n_subdivisions)
3591 for (
unsigned int data_set = 0; data_set < n_data_sets;
3593 out << patch.data(data_set, i1 * d1 + i2 * d2) <<
' ';
3609 for (
unsigned int i3 = 0; i3 < n3; ++i3)
3610 for (
unsigned int i2 = 0; i2 < n2; ++i2)
3611 for (
unsigned int i1 = 0; i1 < n1; ++i1)
3615 compute_node(patch, i1, i2, i3, n_subdivisions);
3617 if (i1 < n_subdivisions)
3621 for (
unsigned int data_set = 0; data_set < n_data_sets;
3624 << patch.data(data_set,
3625 i1 * d1 + i2 * d2 + i3 * d3);
3629 out << compute_node(
3630 patch, i1 + 1, i2, i3, n_subdivisions);
3632 for (
unsigned int data_set = 0; data_set < n_data_sets;
3635 << patch.data(data_set,
3636 (i1 + 1) * d1 + i2 * d2 + i3 * d3);
3640 out <<
'\n' <<
'\n';
3644 if (i2 < n_subdivisions)
3648 for (
unsigned int data_set = 0; data_set < n_data_sets;
3651 << patch.data(data_set,
3652 i1 * d1 + i2 * d2 + i3 * d3);
3656 out << compute_node(
3657 patch, i1, i2 + 1, i3, n_subdivisions);
3659 for (
unsigned int data_set = 0; data_set < n_data_sets;
3662 << patch.data(data_set,
3663 i1 * d1 + (i2 + 1) * d2 + i3 * d3);
3667 out <<
'\n' <<
'\n';
3671 if (i3 < n_subdivisions)
3675 for (
unsigned int data_set = 0; data_set < n_data_sets;
3678 << patch.data(data_set,
3679 i1 * d1 + i2 * d2 + i3 * d3);
3683 out << compute_node(
3684 patch, i1, i2, i3 + 1, n_subdivisions);
3686 for (
unsigned int data_set = 0; data_set < n_data_sets;
3689 << patch.data(data_set,
3690 i1 * d1 + i2 * d2 + (i3 + 1) * d3);
3693 out <<
'\n' <<
'\n';
3708 template <
int dim,
int spacedim>
3712 const std::vector<std::string> & data_names,
3714 std::tuple<
unsigned int,
3723 #ifndef DEAL_II_WITH_MPI 3732 if (patches.size() == 0)
3739 const unsigned int n_data_sets = data_names.size();
3744 out <<
"/* This file was generated by the deal.II library." <<
'\n' 3748 <<
" For a description of the POVRAY format see the POVRAY manual." 3753 out <<
"#include \"colors.inc\" " <<
'\n' 3754 <<
"#include \"textures.inc\" " <<
'\n';
3759 out <<
"#include \"data.inc\" " <<
'\n';
3765 <<
"camera {" <<
'\n' 3766 <<
" location <1,4,-7>" <<
'\n' 3767 <<
" look_at <0,0,0>" <<
'\n' 3768 <<
" angle 30" <<
'\n' 3773 <<
"light_source {" <<
'\n' 3774 <<
" <1,4,-7>" <<
'\n' 3775 <<
" color Grey" <<
'\n' 3778 <<
"light_source {" <<
'\n' 3779 <<
" <0,20,0>" <<
'\n' 3780 <<
" color White" <<
'\n' 3787 double hmin = patches[0].data(0, 0);
3788 double hmax = patches[0].data(0, 0);
3790 for (
const auto &patch : patches)
3792 const unsigned int n_subdivisions = patch.n_subdivisions;
3794 Assert((patch.data.n_rows() == n_data_sets &&
3795 !patch.points_are_available) ||
3796 (patch.data.n_rows() == n_data_sets + spacedim &&
3797 patch.points_are_available),
3799 (n_data_sets + spacedim) :
3801 patch.data.n_rows()));
3802 Assert(patch.data.n_cols() ==
3803 Utilities::fixed_power<dim>(n_subdivisions + 1),
3806 for (
unsigned int i = 0; i < n_subdivisions + 1; ++i)
3807 for (
unsigned int j = 0; j < n_subdivisions + 1; ++j)
3809 const int dl = i * (n_subdivisions + 1) + j;
3810 if (patch.data(0, dl) < hmin)
3811 hmin = patch.data(0, dl);
3812 if (patch.data(0, dl) > hmax)
3813 hmax = patch.data(0, dl);
3817 out <<
"#declare HMIN=" << hmin <<
";" <<
'\n' 3818 <<
"#declare HMAX=" << hmax <<
";" <<
'\n' 3824 out <<
"#declare Tex=texture{" <<
'\n' 3825 <<
" pigment {" <<
'\n' 3826 <<
" gradient y" <<
'\n' 3827 <<
" scale y*(HMAX-HMIN)*" << 0.1 <<
'\n' 3828 <<
" color_map {" <<
'\n' 3829 <<
" [0.00 color Light_Purple] " <<
'\n' 3830 <<
" [0.95 color Light_Purple] " <<
'\n' 3831 <<
" [1.00 color White] " <<
'\n' 3839 out <<
'\n' <<
"mesh {" <<
'\n';
3843 for (
const auto &patch : patches)
3845 const unsigned int n_subdivisions = patch.n_subdivisions;
3846 const unsigned int n = n_subdivisions + 1;
3847 const unsigned int d1 = 1;
3848 const unsigned int d2 = n;
3850 Assert((patch.data.n_rows() == n_data_sets &&
3851 !patch.points_are_available) ||
3852 (patch.data.n_rows() == n_data_sets + spacedim &&
3853 patch.points_are_available),
3855 (n_data_sets + spacedim) :
3857 patch.data.n_rows()));
3858 Assert(patch.data.n_cols() == Utilities::fixed_power<dim>(n),
3862 std::vector<Point<spacedim>> ver(n * n);
3864 for (
unsigned int i2 = 0; i2 < n; ++i2)
3865 for (
unsigned int i1 = 0; i1 < n; ++i1)
3868 ver[i1 * d1 + i2 * d2] =
3869 compute_node(patch, i1, i2, 0, n_subdivisions);
3876 std::vector<Point<3>> nrml;
3886 for (
unsigned int i = 0; i < n; ++i)
3887 for (
unsigned int j = 0; j < n; ++j)
3889 const unsigned int il = (i == 0) ? i : (i - 1);
3890 const unsigned int ir =
3891 (i == n_subdivisions) ? i : (i + 1);
3892 const unsigned int jl = (j == 0) ? j : (j - 1);
3893 const unsigned int jr =
3894 (j == n_subdivisions) ? j : (j + 1);
3897 ver[ir * d1 + j * d2](0) - ver[il * d1 + j * d2](0);
3898 h1(1) = patch.data(0, ir * d1 + j * d2) -
3899 patch.data(0, il * d1 + j * d2);
3901 ver[ir * d1 + j * d2](1) - ver[il * d1 + j * d2](1);
3904 ver[i * d1 + jr * d2](0) - ver[i * d1 + jl * d2](0);
3905 h2(1) = patch.data(0, i * d1 + jr * d2) -
3906 patch.data(0, i * d1 + jl * d2);
3908 ver[i * d1 + jr * d2](1) - ver[i * d1 + jl * d2](1);
3910 nrml[i * d1 + j * d2](0) = h1(1) * h2(2) - h1(2) * h2(1);
3911 nrml[i * d1 + j * d2](1) = h1(2) * h2(0) - h1(0) * h2(2);
3912 nrml[i * d1 + j * d2](2) = h1(0) * h2(1) - h1(1) * h2(0);
3916 std::sqrt(std::pow(nrml[i * d1 + j * d2](0), 2.) +
3917 std::pow(nrml[i * d1 + j * d2](1), 2.) +
3918 std::pow(nrml[i * d1 + j * d2](2), 2.));
3920 if (nrml[i * d1 + j * d2](1) < 0)
3923 for (
unsigned int k = 0; k < 3; ++k)
3924 nrml[i * d1 + j * d2](k) /=
norm;
3929 for (
unsigned int i = 0; i < n_subdivisions; ++i)
3930 for (
unsigned int j = 0; j < n_subdivisions; ++j)
3933 const int dl = i * d1 + j * d2;
3939 out <<
"smooth_triangle {" <<
'\n' 3940 <<
"\t<" << ver[dl](0) <<
"," << patch.data(0, dl)
3941 <<
"," << ver[dl](1) <<
">, <" << nrml[dl](0) <<
", " 3942 << nrml[dl](1) <<
", " << nrml[dl](2) <<
">," <<
'\n';
3943 out <<
" \t<" << ver[dl + d1](0) <<
"," 3944 << patch.data(0, dl + d1) <<
"," << ver[dl + d1](1)
3945 <<
">, <" << nrml[dl + d1](0) <<
", " 3946 << nrml[dl + d1](1) <<
", " << nrml[dl + d1](2)
3948 out <<
"\t<" << ver[dl + d1 + d2](0) <<
"," 3949 << patch.data(0, dl + d1 + d2) <<
"," 3950 << ver[dl + d1 + d2](1) <<
">, <" 3951 << nrml[dl + d1 + d2](0) <<
", " 3952 << nrml[dl + d1 + d2](1) <<
", " 3953 << nrml[dl + d1 + d2](2) <<
">}" <<
'\n';
3956 out <<
"smooth_triangle {" <<
'\n' 3957 <<
"\t<" << ver[dl](0) <<
"," << patch.data(0, dl)
3958 <<
"," << ver[dl](1) <<
">, <" << nrml[dl](0) <<
", " 3959 << nrml[dl](1) <<
", " << nrml[dl](2) <<
">," <<
'\n';
3960 out <<
"\t<" << ver[dl + d1 + d2](0) <<
"," 3961 << patch.data(0, dl + d1 + d2) <<
"," 3962 << ver[dl + d1 + d2](1) <<
">, <" 3963 << nrml[dl + d1 + d2](0) <<
", " 3964 << nrml[dl + d1 + d2](1) <<
", " 3965 << nrml[dl + d1 + d2](2) <<
">," <<
'\n';
3966 out <<
"\t<" << ver[dl + d2](0) <<
"," 3967 << patch.data(0, dl + d2) <<
"," << ver[dl + d2](1)
3968 <<
">, <" << nrml[dl + d2](0) <<
", " 3969 << nrml[dl + d2](1) <<
", " << nrml[dl + d2](2)
3975 out <<
"triangle {" <<
'\n' 3976 <<
"\t<" << ver[dl](0) <<
"," << patch.data(0, dl)
3977 <<
"," << ver[dl](1) <<
">," <<
'\n';
3978 out <<
"\t<" << ver[dl + d1](0) <<
"," 3979 << patch.data(0, dl + d1) <<
"," << ver[dl + d1](1)
3981 out <<
"\t<" << ver[dl + d1 + d2](0) <<
"," 3982 << patch.data(0, dl + d1 + d2) <<
"," 3983 << ver[dl + d1 + d2](1) <<
">}" <<
'\n';
3986 out <<
"triangle {" <<
'\n' 3987 <<
"\t<" << ver[dl](0) <<
"," << patch.data(0, dl)
3988 <<
"," << ver[dl](1) <<
">," <<
'\n';
3989 out <<
"\t<" << ver[dl + d1 + d2](0) <<
"," 3990 << patch.data(0, dl + d1 + d2) <<
"," 3991 << ver[dl + d1 + d2](1) <<
">," <<
'\n';
3992 out <<
"\t<" << ver[dl + d2](0) <<
"," 3993 << patch.data(0, dl + d2) <<
"," << ver[dl + d2](1)
4001 Assert(n_subdivisions == 3,
4004 <<
"bicubic_patch {" <<
'\n' 4005 <<
" type 0" <<
'\n' 4006 <<
" flatness 0" <<
'\n' 4007 <<
" u_steps 0" <<
'\n' 4008 <<
" v_steps 0" <<
'\n';
4009 for (
int i = 0; i < 16; ++i)
4011 out <<
"\t<" << ver[i](0) <<
"," << patch.data(0, i) <<
"," 4012 << ver[i](1) <<
">";
4017 out <<
" texture {Tex}" <<
'\n' <<
"}" <<
'\n';
4024 out <<
" texture {Tex}" <<
'\n' <<
"}" <<
'\n' <<
'\n';
4035 template <
int dim,
int spacedim>
4039 const std::vector<std::string> & ,
4041 std::tuple<
unsigned int,
4053 template <
int spacedim>
4057 const std::vector<std::string> & ,
4059 std::tuple<
unsigned int,
4068 #ifndef DEAL_II_WITH_MPI 4077 if (patches.size() == 0)
4087 std::multiset<EpsCell2d> cells;
4096 double heights[4] = {0, 0, 0, 0};
4100 for (
const auto &patch : patches)
4102 const unsigned int n_subdivisions = patch.n_subdivisions;
4103 const unsigned int n = n_subdivisions + 1;
4104 const unsigned int d1 = 1;
4105 const unsigned int d2 = n;
4107 for (
unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
4108 for (
unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
4111 points[0] = compute_node(patch, i1, i2, 0, n_subdivisions);
4112 points[1] = compute_node(patch, i1 + 1, i2, 0, n_subdivisions);
4113 points[2] = compute_node(patch, i1, i2 + 1, 0, n_subdivisions);
4115 compute_node(patch, i1 + 1, i2 + 1, 0, n_subdivisions);
4121 patch.data.n_rows() == 0,
4124 patch.data.n_rows()));
4126 patch.data.n_rows() != 0 ?
4130 heights[1] = patch.data.n_rows() != 0 ?
4132 (i1 + 1) * d1 + i2 * d2) *
4135 heights[2] = patch.data.n_rows() != 0 ?
4137 i1 * d1 + (i2 + 1) * d2) *
4140 heights[3] = patch.data.n_rows() != 0 ?
4142 (i1 + 1) * d1 + (i2 + 1) * d2) *
4149 for (
unsigned int i = 0; i < 4; ++i)
4150 heights[i] = points[i](2);
4171 cz = -std::cos(flags.
turn_angle * 2 * pi / 360.),
4174 sz = std::sin(flags.
turn_angle * 2 * pi / 360.);
4175 for (
unsigned int vertex = 0; vertex < 4; ++vertex)
4177 const double x = points[vertex](0), y = points[vertex](1),
4178 z = -heights[vertex];
4180 eps_cell.vertices[vertex](0) = -cz * x + sz * y;
4181 eps_cell.vertices[vertex](1) =
4182 -cx * sz * x - cx * cz * y - sx * z;
4206 (points[0] + points[1] + points[2] + points[3]) / 4;
4207 const double center_height =
4208 -(heights[0] + heights[1] + heights[2] + heights[3]) / 4;
4211 eps_cell.depth = -sx * sz * center_point(0) -
4212 sx * cz * center_point(1) + cx * center_height;
4217 patch.data.n_rows() == 0,
4220 patch.data.n_rows()));
4221 const double color_values[4] = {
4222 patch.data.n_rows() != 0 ?
4226 patch.data.n_rows() != 0 ?
4227 patch.data(flags.
color_vector, (i1 + 1) * d1 + i2 * d2) :
4230 patch.data.n_rows() != 0 ?
4231 patch.data(flags.
color_vector, i1 * d1 + (i2 + 1) * d2) :
4234 patch.data.n_rows() != 0 ?
4236 (i1 + 1) * d1 + (i2 + 1) * d2) :
4240 eps_cell.color_value = (color_values[0] + color_values[1] +
4241 color_values[3] + color_values[2]) /
4246 std::min(min_color_value, eps_cell.color_value);
4248 std::max(max_color_value, eps_cell.color_value);
4252 cells.insert(eps_cell);
4258 double x_min = cells.begin()->vertices[0](0);
4259 double x_max = x_min;
4260 double y_min = cells.begin()->vertices[0](1);
4261 double y_max = y_min;
4263 for (
const auto &cell : cells)
4264 for (
const auto &vertex : cell.vertices)
4266 x_min =
std::min(x_min, vertex(0));
4267 x_max =
std::max(x_max, vertex(0));
4268 y_min =
std::min(y_min, vertex(1));
4269 y_max =
std::max(y_max, vertex(1));
4274 const double scale =
4278 const Point<2> offset(x_min, y_min);
4283 out <<
"%!PS-Adobe-2.0 EPSF-1.2" <<
'\n' 4284 <<
"%%Title: deal.II Output" <<
'\n' 4285 <<
"%%Creator: the deal.II library" <<
'\n' 4288 <<
"%%BoundingBox: " 4292 <<
static_cast<unsigned int>((x_max - x_min) * scale + 0.5) <<
' ' 4293 <<
static_cast<unsigned int>((y_max - y_min) * scale + 0.5) <<
'\n';
4302 out <<
"/m {moveto} bind def" <<
'\n' 4303 <<
"/l {lineto} bind def" <<
'\n' 4304 <<
"/s {setrgbcolor} bind def" <<
'\n' 4305 <<
"/sg {setgray} bind def" <<
'\n' 4306 <<
"/lx {lineto closepath stroke} bind def" <<
'\n' 4307 <<
"/lf {lineto closepath fill} bind def" <<
'\n';
4309 out <<
"%%EndProlog" <<
'\n' <<
'\n';
4311 out << flags.
line_width <<
" setlinewidth" <<
'\n';
4319 if (max_color_value == min_color_value)
4320 max_color_value = min_color_value + 1;
4324 for (
const auto &cell : cells)
4337 out << rgb_values.
red <<
" sg ";
4339 out << rgb_values.
red <<
' ' << rgb_values.
green <<
' ' 4340 << rgb_values.
blue <<
" s ";
4345 out << (cell.vertices[0] - offset) * scale <<
" m " 4346 << (cell.vertices[1] - offset) * scale <<
" l " 4347 << (cell.vertices[3] - offset) * scale <<
" l " 4348 << (cell.vertices[2] - offset) * scale <<
" lf" <<
'\n';
4353 << (cell.vertices[0] - offset) * scale <<
" m " 4354 << (cell.vertices[1] - offset) * scale <<
" l " 4355 << (cell.vertices[3] - offset) * scale <<
" l " 4356 << (cell.vertices[2] - offset) * scale <<
" lx" <<
'\n';
4358 out <<
"showpage" <<
'\n';
4367 template <
int dim,
int spacedim>
4371 const std::vector<std::string> & data_names,
4373 std::tuple<
unsigned int,
4389 #ifndef DEAL_II_WITH_MPI 4398 if (patches.size() == 0)
4402 GmvStream gmv_out(out, flags);
4403 const unsigned int n_data_sets = data_names.size();
4406 Assert((patches[0].data.n_rows() == n_data_sets &&
4407 !patches[0].points_are_available) ||
4408 (patches[0].data.n_rows() == n_data_sets + spacedim &&
4409 patches[0].points_are_available),
4411 (n_data_sets + spacedim) :
4413 patches[0].data.n_rows()));
4417 out <<
"gmvinput ascii" <<
'\n' <<
'\n';
4420 unsigned int n_nodes;
4422 compute_sizes<dim, spacedim>(patches, n_nodes,
n_cells);
4437 void (*fun_ptr)(
const std::vector<Patch<dim, spacedim>> &,
4439 &write_gmv_reorder_data_vectors<dim, spacedim>;
4447 out <<
"nodes " << n_nodes <<
'\n';
4448 for (
unsigned int d = 0;
d < spacedim; ++
d)
4450 gmv_out.selected_component =
d;
4456 for (
unsigned int d = spacedim;
d < 3; ++
d)
4458 for (
unsigned int i = 0; i < n_nodes; ++i)
4465 out <<
"cells " << n_cells <<
'\n';
4470 out <<
"variable" <<
'\n';
4474 reorder_task.
join();
4478 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
4480 out << data_names[data_set] <<
" 1" <<
'\n';
4482 data_vectors[data_set].
end(),
4483 std::ostream_iterator<double>(out,
" "));
4484 out <<
'\n' <<
'\n';
4490 out <<
"endvars" <<
'\n';
4493 out <<
"endgmv" <<
'\n';
4504 template <
int dim,
int spacedim>
4508 const std::vector<std::string> & data_names,
4510 std::tuple<
unsigned int,
4524 #ifndef DEAL_II_WITH_MPI 4533 if (patches.size() == 0)
4537 TecplotStream tecplot_out(out, flags);
4539 const unsigned int n_data_sets = data_names.size();
4542 Assert((patches[0].data.n_rows() == n_data_sets &&
4543 !patches[0].points_are_available) ||
4544 (patches[0].data.n_rows() == n_data_sets + spacedim &&
4545 patches[0].points_are_available),
4547 (n_data_sets + spacedim) :
4549 patches[0].data.n_rows()));
4552 unsigned int n_nodes;
4554 compute_sizes<dim, spacedim>(patches, n_nodes,
n_cells);
4560 <<
"# This file was generated by the deal.II library." <<
'\n' 4564 <<
"# For a description of the Tecplot format see the Tecplot documentation." 4569 out <<
"Variables=";
4577 out <<
"\"x\", \"y\"";
4580 out <<
"\"x\", \"y\", \"z\"";
4586 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
4587 out <<
", \"" << data_names[data_set] <<
"\"";
4593 out <<
"t=\"" << flags.
zone_name <<
"\" ";
4596 out <<
"strandid=1, solutiontime=" << flags.
solution_time <<
", ";
4598 out <<
"f=feblock, n=" << n_nodes <<
", e=" << n_cells
4599 <<
", et=" << tecplot_cell_type[dim] <<
'\n';
4618 void (*fun_ptr)(
const std::vector<Patch<dim, spacedim>> &,
4620 &write_gmv_reorder_data_vectors<dim, spacedim>;
4628 for (
unsigned int d = 0;
d < spacedim; ++
d)
4630 tecplot_out.selected_component =
d;
4641 reorder_task.
join();
4644 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
4647 data_vectors[data_set].
end(),
4648 std::ostream_iterator<double>(out,
"\n"));
4666 #ifdef DEAL_II_HAVE_TECPLOT 4673 TecplotMacros(
const unsigned int n_nodes = 0,
4674 const unsigned int n_vars = 0,
4675 const unsigned int n_cells = 0,
4676 const unsigned int n_vert = 0);
4679 nd(
const unsigned int i,
const unsigned int j);
4681 cd(
const unsigned int i,
const unsigned int j);
4682 std::vector<float> nodalData;
4683 std::vector<int> connData;
4686 unsigned int n_nodes;
4687 unsigned int n_vars;
4689 unsigned int n_vert;
4693 inline TecplotMacros::TecplotMacros(
const unsigned int n_nodes,
4694 const unsigned int n_vars,
4696 const unsigned int n_vert)
4702 nodalData.resize(n_nodes * n_vars);
4703 connData.resize(n_cells * n_vert);
4708 inline TecplotMacros::~TecplotMacros()
4714 TecplotMacros::nd(
const unsigned int i,
const unsigned int j)
4716 return nodalData[i * n_nodes + j];
4722 TecplotMacros::cd(
const unsigned int i,
const unsigned int j)
4724 return connData[i + j * n_vert];
4735 template <
int dim,
int spacedim>
4739 const std::vector<std::string> & data_names,
4741 std::tuple<
unsigned int,
4745 & nonscalar_data_ranges,
4754 #ifndef DEAL_II_HAVE_TECPLOT 4757 write_tecplot(patches, data_names, nonscalar_data_ranges, flags, out);
4765 write_tecplot(patches, data_names, nonscalar_data_ranges, flags, out);
4772 char *file_name = (
char *)flags.tecplot_binary_file_name;
4774 if (file_name ==
nullptr)
4779 ExcMessage(
"Specify the name of the tecplot_binary" 4780 " file through the TecplotFlags interface."));
4781 write_tecplot(patches, data_names, nonscalar_data_ranges, flags, out);
4788 # ifndef DEAL_II_WITH_MPI 4797 if (patches.size() == 0)
4801 const unsigned int n_data_sets = data_names.size();
4804 Assert((patches[0].data.n_rows() == n_data_sets &&
4805 !patches[0].points_are_available) ||
4806 (patches[0].data.n_rows() == n_data_sets + spacedim &&
4807 patches[0].points_are_available),
4809 (n_data_sets + spacedim) :
4811 patches[0].data.n_rows()));
4814 unsigned int n_nodes;
4816 compute_sizes<dim, spacedim>(patches, n_nodes,
n_cells);
4818 const unsigned int vars_per_node = (spacedim + n_data_sets),
4821 TecplotMacros tm(n_nodes, vars_per_node, n_cells, nodes_per_cell);
4823 int is_double = 0, tec_debug = 0, cell_type = tecplot_binary_cell_type[dim];
4825 std::string tec_var_names;
4829 tec_var_names =
"x y";
4832 tec_var_names =
"x y z";
4838 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
4840 tec_var_names +=
" ";
4841 tec_var_names += data_names[data_set];
4857 void (*fun_ptr)(
const std::vector<Patch<dim, spacedim>> &,
4859 &write_gmv_reorder_data_vectors<dim, spacedim>;
4865 for (
unsigned int d = 1;
d <= spacedim; ++
d)
4867 unsigned int entry = 0;
4869 for (
const auto &patch : patches)
4871 const unsigned int n_subdivisions = patch.n_subdivisions;
4877 for (
unsigned int j = 0; j < n_subdivisions + 1; ++j)
4878 for (
unsigned int i = 0; i < n_subdivisions + 1; ++i)
4880 const double x_frac = i * 1. / n_subdivisions,
4881 y_frac = j * 1. / n_subdivisions;
4883 tm.nd((
d - 1), entry) =
static_cast<float>(
4884 (((patch.vertices[1](
d - 1) * x_frac) +
4885 (patch.vertices[0](
d - 1) * (1 - x_frac))) *
4887 ((patch.vertices[3](
d - 1) * x_frac) +
4888 (patch.vertices[2](
d - 1) * (1 - x_frac))) *
4897 for (
unsigned int j = 0; j < n_subdivisions + 1; ++j)
4898 for (
unsigned int k = 0; k < n_subdivisions + 1; ++k)
4899 for (
unsigned int i = 0; i < n_subdivisions + 1; ++i)
4901 const double x_frac = i * 1. / n_subdivisions,
4902 y_frac = k * 1. / n_subdivisions,
4903 z_frac = j * 1. / n_subdivisions;
4906 tm.nd((
d - 1), entry) =
static_cast<float>(
4907 ((((patch.vertices[1](
d - 1) * x_frac) +
4908 (patch.vertices[0](
d - 1) * (1 - x_frac))) *
4910 ((patch.vertices[3](
d - 1) * x_frac) +
4911 (patch.vertices[2](
d - 1) * (1 - x_frac))) *
4914 (((patch.vertices[5](
d - 1) * x_frac) +
4915 (patch.vertices[4](
d - 1) * (1 - x_frac))) *
4917 ((patch.vertices[7](
d - 1) * x_frac) +
4918 (patch.vertices[6](
d - 1) * (1 - x_frac))) *
4936 reorder_task.
join();
4939 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
4940 for (
unsigned int entry = 0; entry < data_vectors[data_set].
size();
4942 tm.nd((spacedim + data_set), entry) =
4943 static_cast<float>(data_vectors[data_set][entry]);
4949 unsigned int first_vertex_of_patch = 0;
4950 unsigned int elem = 0;
4952 for (
const auto &patch : patches)
4954 const unsigned int n_subdivisions = patch.n_subdivisions;
4955 const unsigned int n = n_subdivisions + 1;
4956 const unsigned int d1 = 1;
4957 const unsigned int d2 = n;
4958 const unsigned int d3 = n * n;
4964 for (
unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
4965 for (
unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
4968 first_vertex_of_patch + (i1)*d1 + (i2)*d2 + 1;
4970 first_vertex_of_patch + (i1 + 1) * d1 + (i2)*d2 + 1;
4971 tm.cd(2, elem) = first_vertex_of_patch + (i1 + 1) * d1 +
4974 first_vertex_of_patch + (i1)*d1 + (i2 + 1) * d2 + 1;
4983 for (
unsigned int i3 = 0; i3 < n_subdivisions; ++i3)
4984 for (
unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
4985 for (
unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
4990 tm.cd(0, elem) = first_vertex_of_patch + (i1)*d1 +
4991 (i2)*d2 + (i3)*d3 + 1;
4992 tm.cd(1, elem) = first_vertex_of_patch + (i1 + 1) * d1 +
4993 (i2)*d2 + (i3)*d3 + 1;
4994 tm.cd(2, elem) = first_vertex_of_patch + (i1 + 1) * d1 +
4995 (i2 + 1) * d2 + (i3)*d3 + 1;
4996 tm.cd(3, elem) = first_vertex_of_patch + (i1)*d1 +
4997 (i2 + 1) * d2 + (i3)*d3 + 1;
4998 tm.cd(4, elem) = first_vertex_of_patch + (i1)*d1 +
4999 (i2)*d2 + (i3 + 1) * d3 + 1;
5000 tm.cd(5, elem) = first_vertex_of_patch + (i1 + 1) * d1 +
5001 (i2)*d2 + (i3 + 1) * d3 + 1;
5002 tm.cd(6, elem) = first_vertex_of_patch + (i1 + 1) * d1 +
5003 (i2 + 1) * d2 + (i3 + 1) * d3 + 1;
5004 tm.cd(7, elem) = first_vertex_of_patch + (i1)*d1 +
5005 (i2 + 1) * d2 + (i3 + 1) * d3 + 1;
5018 first_vertex_of_patch += Utilities::fixed_power<dim>(n);
5023 int ierr = 0, num_nodes =
static_cast<int>(n_nodes),
5024 num_cells = static_cast<int>(n_cells);
5026 char dot[2] = {
'.', 0};
5030 char *var_names =
const_cast<char *
>(tec_var_names.c_str());
5031 ierr = TECINI(
nullptr, var_names, file_name, dot, &tec_debug, &is_double);
5035 char FEBLOCK[] = {
'F',
'E',
'B',
'L',
'O',
'C',
'K', 0};
5037 TECZNE(
nullptr, &num_nodes, &num_cells, &cell_type, FEBLOCK,
nullptr);
5041 int total = (vars_per_node * num_nodes);
5043 ierr = TECDAT(&total, tm.nodalData.data(), &is_double);
5047 ierr = TECNOD(tm.connData.data());
5060 template <
int dim,
int spacedim>
5064 const std::vector<std::string> & data_names,
5066 std::tuple<
unsigned int,
5070 & nonscalar_data_ranges,
5076 #ifndef DEAL_II_WITH_MPI 5085 if (patches.size() == 0)
5089 VtkStream vtk_out(out, flags);
5091 const unsigned int n_data_sets = data_names.size();
5093 if (patches[0].points_are_available)
5105 out <<
"# vtk DataFile Version 3.0" <<
'\n' 5106 <<
"#This file was generated by the deal.II library";
5114 out <<
'\n' <<
"ASCII" <<
'\n';
5116 out <<
"DATASET UNSTRUCTURED_GRID\n" <<
'\n';
5123 const unsigned int n_metadata =
5128 out <<
"FIELD FieldData " << n_metadata <<
"\n";
5132 out <<
"CYCLE 1 1 int\n" << flags.
cycle <<
"\n";
5136 out <<
"TIME 1 1 double\n" << flags.
time <<
"\n";
5142 unsigned int n_nodes,
n_cells, n_points_and_n_cells;
5143 compute_sizes(patches,
5147 n_points_and_n_cells);
5163 void (*fun_ptr)(
const std::vector<Patch<dim, spacedim>> &,
5165 &write_gmv_reorder_data_vectors<dim, spacedim>;
5173 out <<
"POINTS " << n_nodes <<
" double" <<
'\n';
5178 out <<
"CELLS " << n_cells <<
' ' << n_points_and_n_cells <<
'\n';
5186 out <<
"CELL_TYPES " << n_cells <<
'\n';
5190 for (
const auto &patch : patches)
5192 const auto vtk_cell_id =
5195 for (
unsigned int i = 0; i < vtk_cell_id[1]; ++i)
5196 out <<
' ' << vtk_cell_id[0];
5205 reorder_task.
join();
5210 out <<
"POINT_DATA " << n_nodes <<
'\n';
5214 std::vector<bool> data_set_written(n_data_sets,
false);
5215 for (
const auto &nonscalar_data_range : nonscalar_data_ranges)
5222 std::get<0>(nonscalar_data_range),
5224 std::get<0>(nonscalar_data_range)));
5225 AssertThrow(std::get<1>(nonscalar_data_range) < n_data_sets,
5229 AssertThrow(std::get<1>(nonscalar_data_range) + 1 -
5230 std::get<0>(nonscalar_data_range) <=
5233 "Can't declare a vector with more than 3 components " 5237 for (
unsigned int i = std::get<0>(nonscalar_data_range);
5238 i <= std::get<1>(nonscalar_data_range);
5240 data_set_written[i] =
true;
5246 if (!std::get<2>(nonscalar_data_range).empty())
5247 out << std::get<2>(nonscalar_data_range);
5250 for (
unsigned int i = std::get<0>(nonscalar_data_range);
5251 i < std::get<1>(nonscalar_data_range);
5253 out << data_names[i] <<
"__";
5254 out << data_names[std::get<1>(nonscalar_data_range)];
5257 out <<
" double" <<
'\n';
5260 for (
unsigned int n = 0; n < n_nodes; ++n)
5262 switch (std::get<1>(nonscalar_data_range) -
5263 std::get<0>(nonscalar_data_range))
5266 out << data_vectors(std::get<0>(nonscalar_data_range), n)
5271 out << data_vectors(std::get<0>(nonscalar_data_range), n)
5273 << data_vectors(std::get<0>(nonscalar_data_range) + 1, n)
5277 out << data_vectors(std::get<0>(nonscalar_data_range), n)
5279 << data_vectors(std::get<0>(nonscalar_data_range) + 1, n)
5281 << data_vectors(std::get<0>(nonscalar_data_range) + 2, n)
5294 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
5295 if (data_set_written[data_set] ==
false)
5297 out <<
"SCALARS " << data_names[data_set] <<
" double 1" <<
'\n' 5298 <<
"LOOKUP_TABLE default" <<
'\n';
5300 data_vectors[data_set].
end(),
5301 std::ostream_iterator<double>(out,
" "));
5317 out <<
"<?xml version=\"1.0\" ?> \n";
5319 out <<
"# vtk DataFile Version 3.0" <<
'\n' 5320 <<
"#This file was generated by the deal.II library";
5329 out <<
"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\"";
5330 #ifdef DEAL_II_WITH_ZLIB 5331 out <<
" compressor=\"vtkZLibDataCompressor\"";
5333 #ifdef DEAL_II_WORDS_BIGENDIAN 5334 out <<
" byte_order=\"BigEndian\"";
5336 out <<
" byte_order=\"LittleEndian\"";
5340 out <<
"<UnstructuredGrid>";
5350 out <<
" </UnstructuredGrid>\n";
5351 out <<
"</VTKFile>\n";
5356 template <
int dim,
int spacedim>
5360 const std::vector<std::string> & data_names,
5362 std::tuple<
unsigned int,
5366 & nonscalar_data_ranges,
5371 write_vtu_main(patches, data_names, nonscalar_data_ranges, flags, out);
5378 template <
int dim,
int spacedim>
5382 const std::vector<std::string> & data_names,
5384 std::tuple<
unsigned int,
5388 & nonscalar_data_ranges,
5394 #ifndef DEAL_II_WITH_MPI 5403 if (patches.size() == 0)
5408 out <<
"<Piece NumberOfPoints=\"0\" NumberOfCells=\"0\" >\n" 5410 <<
"<DataArray type=\"UInt8\" Name=\"types\"></DataArray>\n" 5412 <<
" <PointData Scalars=\"scalars\">\n";
5413 std::vector<bool> data_set_written(data_names.size(),
false);
5414 for (
const auto &nonscalar_data_range : nonscalar_data_ranges)
5417 for (
unsigned int i = std::get<0>(nonscalar_data_range);
5418 i <= std::get<1>(nonscalar_data_range);
5420 data_set_written[i] =
true;
5424 out <<
" <DataArray type=\"Float32\" Name=\"";
5426 if (!std::get<2>(nonscalar_data_range).empty())
5427 out << std::get<2>(nonscalar_data_range);
5430 for (
unsigned int i = std::get<0>(nonscalar_data_range);
5431 i < std::get<1>(nonscalar_data_range);
5433 out << data_names[i] <<
"__";
5434 out << data_names[std::get<1>(nonscalar_data_range)];
5437 out <<
"\" NumberOfComponents=\"3\"></DataArray>\n";
5440 for (
unsigned int data_set = 0; data_set < data_names.size();
5442 if (data_set_written[data_set] ==
false)
5444 out <<
" <DataArray type=\"Float32\" Name=\"" 5445 << data_names[data_set] <<
"\"></DataArray>\n";
5448 out <<
" </PointData>\n";
5449 out <<
"</Piece>\n";
5463 const unsigned int n_metadata =
5467 out <<
"<FieldData>\n";
5472 <<
"<DataArray type=\"Float32\" Name=\"CYCLE\" NumberOfTuples=\"1\" format=\"ascii\">" 5473 << flags.
cycle <<
"</DataArray>\n";
5478 <<
"<DataArray type=\"Float32\" Name=\"TIME\" NumberOfTuples=\"1\" format=\"ascii\">" 5479 << flags.
time <<
"</DataArray>\n";
5483 out <<
"</FieldData>\n";
5487 VtuStream vtu_out(out, flags);
5489 const unsigned int n_data_sets = data_names.size();
5492 if (patches[0].points_are_available)
5501 #ifdef DEAL_II_WITH_ZLIB
5502 const char *ascii_or_binary =
"binary";
5504 const char * ascii_or_binary =
"ascii";
5509 unsigned int n_nodes,
n_cells, n_points_and_n_cells;
5510 compute_sizes(patches,
5514 n_points_and_n_cells);
5530 void (*fun_ptr)(
const std::vector<Patch<dim, spacedim>> &,
5532 &write_gmv_reorder_data_vectors<dim, spacedim, float>;
5541 out <<
"<Piece NumberOfPoints=\"" << n_nodes <<
"\" NumberOfCells=\"" 5542 << n_cells <<
"\" >\n";
5543 out <<
" <Points>\n";
5544 out <<
" <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"" 5545 << ascii_or_binary <<
"\">\n";
5547 out <<
" </DataArray>\n";
5548 out <<
" </Points>\n\n";
5551 out <<
" <Cells>\n";
5552 out <<
" <DataArray type=\"Int32\" Name=\"connectivity\" format=\"" 5553 << ascii_or_binary <<
"\">\n";
5558 out <<
" </DataArray>\n";
5562 out <<
" <DataArray type=\"Int32\" Name=\"offsets\" format=\"" 5563 << ascii_or_binary <<
"\">\n";
5565 std::vector<int32_t> offsets;
5566 offsets.reserve(n_cells);
5570 #ifdef DEAL_II_WITH_ZLIB 5571 std::vector<uint8_t> cell_types;
5573 std::vector<unsigned int> cell_types;
5575 cell_types.reserve(n_cells);
5577 unsigned int first_vertex_of_patch = 0;
5579 for (
const auto &patch : patches)
5581 const auto vtk_cell_id =
5584 for (
unsigned int i = 0; i < vtk_cell_id[1]; ++i)
5586 cell_types.push_back(vtk_cell_id[0]);
5587 first_vertex_of_patch += vtk_cell_id[2];
5588 offsets.push_back(first_vertex_of_patch);
5594 out <<
" </DataArray>\n";
5598 out <<
" <DataArray type=\"UInt8\" Name=\"types\" format=\"" 5599 << ascii_or_binary <<
"\">\n";
5602 vtu_out << cell_types;
5604 out <<
" </DataArray>\n";
5605 out <<
" </Cells>\n";
5613 reorder_task.
join();
5618 out <<
" <PointData Scalars=\"scalars\">\n";
5622 std::vector<bool> data_set_written(n_data_sets,
false);
5623 for (
const auto &range : nonscalar_data_ranges)
5625 const auto first_component = std::get<0>(range);
5626 const auto last_component = std::get<1>(range);
5627 const auto &name = std::get<2>(range);
5628 const bool is_tensor =
5629 (std::get<3>(range) ==
5631 const unsigned int n_components = (is_tensor ? 9 : 3);
5638 AssertThrow((last_component + 1 - first_component <= 9),
5640 "Can't declare a tensor with more than 9 components " 5645 AssertThrow((last_component + 1 - first_component <= 3),
5647 "Can't declare a vector with more than 3 components " 5652 for (
unsigned int i = first_component; i <= last_component; ++i)
5653 data_set_written[i] =
true;
5657 out <<
" <DataArray type=\"Float32\" Name=\"";
5663 for (
unsigned int i = first_component; i < last_component; ++i)
5664 out << data_names[i] <<
"__";
5665 out << data_names[last_component];
5668 out <<
"\" NumberOfComponents=\"" << n_components <<
"\" format=\"" 5669 << ascii_or_binary <<
"\">\n";
5672 std::vector<float> data;
5673 data.reserve(n_nodes * n_components);
5675 for (
unsigned int n = 0; n < n_nodes; ++n)
5679 switch (last_component - first_component)
5682 data.push_back(data_vectors(first_component, n));
5688 data.push_back(data_vectors(first_component, n));
5689 data.push_back(data_vectors(first_component + 1, n));
5694 data.push_back(data_vectors(first_component, n));
5695 data.push_back(data_vectors(first_component + 1, n));
5696 data.push_back(data_vectors(first_component + 2, n));
5709 const unsigned int size = last_component - first_component + 1;
5713 vtk_data[0][0] = data_vectors(first_component, n);
5718 for (
unsigned int c = 0; c < size; ++c)
5722 vtk_data[ind[0]][ind[1]] =
5723 data_vectors(first_component + c, n);
5729 for (
unsigned int c = 0; c < size; ++c)
5733 vtk_data[ind[0]][ind[1]] =
5734 data_vectors(first_component + c, n);
5745 for (
unsigned int i = 0; i < 3; ++i)
5746 for (
unsigned int j = 0; j < 3; ++j)
5747 data.push_back(vtk_data[i][j]);
5752 out <<
" </DataArray>\n";
5757 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
5758 if (data_set_written[data_set] ==
false)
5760 out <<
" <DataArray type=\"Float32\" Name=\"" 5761 << data_names[data_set] <<
"\" format=\"" << ascii_or_binary
5764 std::vector<float> data(data_vectors[data_set].
begin(),
5765 data_vectors[data_set].
end());
5767 out <<
" </DataArray>\n";
5770 out <<
" </PointData>\n";
5773 out <<
" </Piece>\n";
5787 const std::vector<std::string> &piece_names,
5788 const std::vector<std::string> &data_names,
5790 std::tuple<
unsigned int,
5794 &nonscalar_data_ranges)
5798 const unsigned int n_data_sets = data_names.size();
5800 out <<
"<?xml version=\"1.0\"?>\n";
5803 out <<
"#This file was generated by the deal.II library" 5808 <<
"<VTKFile type=\"PUnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
5809 out <<
" <PUnstructuredGrid GhostLevel=\"0\">\n";
5810 out <<
" <PPointData Scalars=\"scalars\">\n";
5813 std::vector<bool> data_set_written(n_data_sets,
false);
5814 for (
const auto &nonscalar_data_range : nonscalar_data_ranges)
5816 const auto first_component = std::get<0>(nonscalar_data_range);
5817 const auto last_component = std::get<1>(nonscalar_data_range);
5818 const bool is_tensor =
5819 (std::get<3>(nonscalar_data_range) ==
5821 const unsigned int n_components = (is_tensor ? 9 : 3);
5828 AssertThrow((last_component + 1 - first_component <= 9),
5830 "Can't declare a tensor with more than 9 components " 5835 Assert((last_component + 1 - first_component <= 3),
5837 "Can't declare a vector with more than 3 components " 5842 for (
unsigned int i = std::get<0>(nonscalar_data_range);
5843 i <= std::get<1>(nonscalar_data_range);
5845 data_set_written[i] =
true;
5849 out <<
" <PDataArray type=\"Float32\" Name=\"";
5851 if (!std::get<2>(nonscalar_data_range).empty())
5852 out << std::get<2>(nonscalar_data_range);
5855 for (
unsigned int i = std::get<0>(nonscalar_data_range);
5856 i < std::get<1>(nonscalar_data_range);
5858 out << data_names[i] <<
"__";
5859 out << data_names[std::get<1>(nonscalar_data_range)];
5862 out <<
"\" NumberOfComponents=\"" << n_components
5863 <<
"\" format=\"ascii\"/>\n";
5866 for (
unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
5867 if (data_set_written[data_set] ==
false)
5869 out <<
" <PDataArray type=\"Float32\" Name=\"" 5870 << data_names[data_set] <<
"\" format=\"ascii\"/>\n";
5873 out <<
" </PPointData>\n";
5875 out <<
" <PPoints>\n";
5876 out <<
" <PDataArray type=\"Float32\" NumberOfComponents=\"3\"/>\n";
5877 out <<
" </PPoints>\n";
5879 for (
const auto &piece_name : piece_names)
5880 out <<
" <Piece Source=\"" << piece_name <<
"\"/>\n";
5882 out <<
" </PUnstructuredGrid>\n";
5883 out <<
"</VTKFile>\n";
5896 const std::vector<std::pair<double, std::string>> ×_and_names)
5900 out <<
"<?xml version=\"1.0\"?>\n";
5903 out <<
"#This file was generated by the deal.II library" 5908 <<
"<VTKFile type=\"Collection\" version=\"0.1\" ByteOrder=\"LittleEndian\">\n";
5909 out <<
" <Collection>\n";
5911 std::streamsize ss = out.precision();
5914 for (
const auto &time_and_name : times_and_names)
5915 out <<
" <DataSet timestep=\"" << time_and_name.first
5916 <<
"\" group=\"\" part=\"0\" file=\"" << time_and_name.second
5919 out <<
" </Collection>\n";
5920 out <<
"</VTKFile>\n";
5932 const std::vector<std::string> &piece_names)
5934 out <<
"!NBLOCKS " << piece_names.size() <<
'\n';
5935 for (
const auto &piece_name : piece_names)
5936 out << piece_name <<
'\n';
5945 const std::vector<std::vector<std::string>> &piece_names)
5949 if (piece_names.size() == 0)
5952 const double nblocks = piece_names[0].size();
5954 ExcMessage(
"piece_names should be a vector of nonempty vectors."));
5956 out <<
"!NBLOCKS " << nblocks <<
'\n';
5957 for (
const auto &domain : piece_names)
5959 Assert(domain.size() == nblocks,
5961 "piece_names should be a vector of equal sized vectors."));
5962 for (
const auto &subdomain : domain)
5963 out << subdomain <<
'\n';
5974 const std::vector<std::pair<
double, std::vector<std::string>>>
5975 ×_and_piece_names)
5979 if (times_and_piece_names.size() == 0)
5982 const double nblocks = times_and_piece_names[0].second.size();
5986 "time_and_piece_names should contain nonempty vectors of filenames for every timestep."));
5988 for (
const auto &domain : times_and_piece_names)
5989 out <<
"!TIME " << domain.first <<
'\n';
5991 out <<
"!NBLOCKS " << nblocks <<
'\n';
5992 for (
const auto &domain : times_and_piece_names)
5994 Assert(domain.second.size() == nblocks,
5996 "piece_names should be a vector of equal sized vectors."));
5997 for (
const auto &subdomain : domain.second)
5998 out << subdomain <<
'\n';
6006 template <
int dim,
int spacedim>
6010 const std::vector<std::string> &,
6012 std::tuple<
unsigned int,
6022 template <
int spacedim>
6026 const std::vector<std::string> & ,
6028 std::tuple<
unsigned int,
6036 const unsigned int height = flags.
height;
6037 unsigned int width = flags.
width;
6040 unsigned int margin_in_percent = 0;
6042 margin_in_percent = 5;
6046 double x_dimension, y_dimension, z_dimension;
6048 const auto &first_patch = patches[0];
6050 unsigned int n_subdivisions = first_patch.n_subdivisions;
6051 unsigned int n = n_subdivisions + 1;
6052 const unsigned int d1 = 1;
6053 const unsigned int d2 = n;
6056 std::array<Point<spacedim>, 4> projected_points;
6059 std::array<Point<2>, 4> projection_decompositions;
6061 projected_point = compute_node(first_patch, 0, 0, 0, n_subdivisions);
6063 if (first_patch.data.n_rows() != 0)
6068 double x_min = projected_point[0];
6069 double x_max = x_min;
6070 double y_min = projected_point[1];
6071 double y_max = y_min;
6072 double z_min = first_patch.data.n_rows() != 0 ?
6075 double z_max = z_min;
6078 for (
const auto &patch : patches)
6080 n_subdivisions = patch.n_subdivisions;
6081 n = n_subdivisions + 1;
6083 for (
unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
6085 for (
unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
6087 projected_points[0] =
6088 compute_node(patch, i1, i2, 0, n_subdivisions);
6089 projected_points[1] =
6090 compute_node(patch, i1 + 1, i2, 0, n_subdivisions);
6091 projected_points[2] =
6092 compute_node(patch, i1, i2 + 1, 0, n_subdivisions);
6093 projected_points[3] =
6094 compute_node(patch, i1 + 1, i2 + 1, 0, n_subdivisions);
6096 x_min =
std::min(x_min, projected_points[0][0]);
6097 x_min =
std::min(x_min, projected_points[1][0]);
6098 x_min =
std::min(x_min, projected_points[2][0]);
6099 x_min =
std::min(x_min, projected_points[3][0]);
6101 x_max =
std::max(x_max, projected_points[0][0]);
6102 x_max =
std::max(x_max, projected_points[1][0]);
6103 x_max =
std::max(x_max, projected_points[2][0]);
6104 x_max =
std::max(x_max, projected_points[3][0]);
6106 y_min =
std::min(y_min, projected_points[0][1]);
6107 y_min =
std::min(y_min, projected_points[1][1]);
6108 y_min =
std::min(y_min, projected_points[2][1]);
6109 y_min =
std::min(y_min, projected_points[3][1]);
6111 y_max =
std::max(y_max, projected_points[0][1]);
6112 y_max =
std::max(y_max, projected_points[1][1]);
6113 y_max =
std::max(y_max, projected_points[2][1]);
6114 y_max =
std::max(y_max, projected_points[3][1]);
6117 patch.data.n_rows() == 0,
6120 patch.data.n_rows()));
6122 z_min = std::min<double>(z_min,
6124 i1 * d1 + i2 * d2));
6125 z_min = std::min<double>(z_min,
6127 (i1 + 1) * d1 + i2 * d2));
6128 z_min = std::min<double>(z_min,
6130 i1 * d1 + (i2 + 1) * d2));
6132 std::min<double>(z_min,
6134 (i1 + 1) * d1 + (i2 + 1) * d2));
6136 z_max = std::max<double>(z_max,
6138 i1 * d1 + i2 * d2));
6139 z_max = std::max<double>(z_max,
6141 (i1 + 1) * d1 + i2 * d2));
6142 z_max = std::max<double>(z_max,
6144 i1 * d1 + (i2 + 1) * d2));
6146 std::max<double>(z_max,
6148 (i1 + 1) * d1 + (i2 + 1) * d2));
6153 x_dimension = x_max - x_min;
6154 y_dimension = y_max - y_min;
6155 z_dimension = z_max - z_min;
6162 float camera_focus = 0;
6165 camera_position[0] = 0.;
6166 camera_position[1] = 0.;
6167 camera_position[2] = z_min + 2. * z_dimension;
6169 camera_direction[0] = 0.;
6170 camera_direction[1] = 0.;
6171 camera_direction[2] = -1.;
6173 camera_horizontal[0] = 1.;
6174 camera_horizontal[1] = 0.;
6175 camera_horizontal[2] = 0.;
6177 camera_focus = .5 * z_dimension;
6183 const float angle_factor = 3.14159265f / 180.f;
6186 camera_position_temp[1] =
6187 std::cos(angle_factor * flags.
polar_angle) * camera_position[1] -
6188 std::sin(angle_factor * flags.
polar_angle) * camera_position[2];
6189 camera_position_temp[2] =
6190 std::sin(angle_factor * flags.
polar_angle) * camera_position[1] +
6191 std::cos(angle_factor * flags.
polar_angle) * camera_position[2];
6193 camera_direction_temp[1] =
6194 std::cos(angle_factor * flags.
polar_angle) * camera_direction[1] -
6195 std::sin(angle_factor * flags.
polar_angle) * camera_direction[2];
6196 camera_direction_temp[2] =
6197 std::sin(angle_factor * flags.
polar_angle) * camera_direction[1] +
6198 std::cos(angle_factor * flags.
polar_angle) * camera_direction[2];
6200 camera_horizontal_temp[1] =
6201 std::cos(angle_factor * flags.
polar_angle) * camera_horizontal[1] -
6202 std::sin(angle_factor * flags.
polar_angle) * camera_horizontal[2];
6203 camera_horizontal_temp[2] =
6204 std::sin(angle_factor * flags.
polar_angle) * camera_horizontal[1] +
6205 std::cos(angle_factor * flags.
polar_angle) * camera_horizontal[2];
6207 camera_position[1] = camera_position_temp[1];
6208 camera_position[2] = camera_position_temp[2];
6210 camera_direction[1] = camera_direction_temp[1];
6211 camera_direction[2] = camera_direction_temp[2];
6213 camera_horizontal[1] = camera_horizontal_temp[1];
6214 camera_horizontal[2] = camera_horizontal_temp[2];
6217 camera_position_temp[0] =
6218 std::cos(angle_factor * flags.
azimuth_angle) * camera_position[0] -
6219 std::sin(angle_factor * flags.
azimuth_angle) * camera_position[1];
6220 camera_position_temp[1] =
6221 std::sin(angle_factor * flags.
azimuth_angle) * camera_position[0] +
6222 std::cos(angle_factor * flags.
azimuth_angle) * camera_position[1];
6224 camera_direction_temp[0] =
6225 std::cos(angle_factor * flags.
azimuth_angle) * camera_direction[0] -
6226 std::sin(angle_factor * flags.
azimuth_angle) * camera_direction[1];
6227 camera_direction_temp[1] =
6228 std::sin(angle_factor * flags.
azimuth_angle) * camera_direction[0] +
6229 std::cos(angle_factor * flags.
azimuth_angle) * camera_direction[1];
6231 camera_horizontal_temp[0] =
6232 std::cos(angle_factor * flags.
azimuth_angle) * camera_horizontal[0] -
6233 std::sin(angle_factor * flags.
azimuth_angle) * camera_horizontal[1];
6234 camera_horizontal_temp[1] =
6235 std::sin(angle_factor * flags.
azimuth_angle) * camera_horizontal[0] +
6236 std::cos(angle_factor * flags.
azimuth_angle) * camera_horizontal[1];
6238 camera_position[0] = camera_position_temp[0];
6239 camera_position[1] = camera_position_temp[1];
6241 camera_direction[0] = camera_direction_temp[0];
6242 camera_direction[1] = camera_direction_temp[1];
6244 camera_horizontal[0] = camera_horizontal_temp[0];
6245 camera_horizontal[1] = camera_horizontal_temp[1];
6248 camera_position[0] = x_min + .5 * x_dimension;
6249 camera_position[1] = y_min + .5 * y_dimension;
6251 camera_position[0] += (z_min + 2. * z_dimension) *
6254 camera_position[1] -= (z_min + 2. * z_dimension) *
6260 double x_min_perspective, y_min_perspective;
6261 double x_max_perspective, y_max_perspective;
6262 double x_dimension_perspective, y_dimension_perspective;
6264 n_subdivisions = first_patch.n_subdivisions;
6265 n = n_subdivisions + 1;
6269 projected_point = compute_node(first_patch, 0, 0, 0, n_subdivisions);
6271 if (first_patch.data.n_rows() != 0)
6276 point[0] = projected_point[0];
6277 point[1] = projected_point[1];
6278 point[2] = first_patch.data.n_rows() != 0 ?
6282 projection_decomposition = svg_project_point(point,
6288 x_min_perspective = projection_decomposition[0];
6289 x_max_perspective = projection_decomposition[0];
6290 y_min_perspective = projection_decomposition[1];
6291 y_max_perspective = projection_decomposition[1];
6294 for (
const auto &patch : patches)
6296 n_subdivisions = patch.n_subdivisions;
6297 for (
unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
6299 for (
unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
6302 {compute_node(patch, i1, i2, 0, n_subdivisions),
6303 compute_node(patch, i1 + 1, i2, 0, n_subdivisions),
6304 compute_node(patch, i1, i2 + 1, 0, n_subdivisions),
6305 compute_node(patch, i1 + 1, i2 + 1, 0, n_subdivisions)}};
6308 patch.data.n_rows() == 0,
6311 patch.data.n_rows()));
6313 const std::array<Point<3>, 4>
vertices = {
6316 patch.data.n_rows() != 0 ?
6317 patch.data(0, i1 * d1 + i2 * d2) :
6321 patch.data.n_rows() != 0 ?
6322 patch.data(0, (i1 + 1) * d1 + i2 * d2) :
6326 patch.data.n_rows() != 0 ?
6327 patch.data(0, i1 * d1 + (i2 + 1) * d2) :
6331 patch.data.n_rows() != 0 ?
6332 patch.data(0, (i1 + 1) * d1 + (i2 + 1) * d2) :
6335 projection_decompositions = {
6336 {svg_project_point(vertices[0],
6341 svg_project_point(vertices[1],
6346 svg_project_point(vertices[2],
6351 svg_project_point(vertices[3],
6359 static_cast<double>(
6360 projection_decompositions[0][0]));
6363 static_cast<double>(
6364 projection_decompositions[1][0]));
6367 static_cast<double>(
6368 projection_decompositions[2][0]));
6371 static_cast<double>(
6372 projection_decompositions[3][0]));
6376 static_cast<double>(
6377 projection_decompositions[0][0]));
6380 static_cast<double>(
6381 projection_decompositions[1][0]));
6384 static_cast<double>(
6385 projection_decompositions[2][0]));
6388 static_cast<double>(
6389 projection_decompositions[3][0]));
6393 static_cast<double>(
6394 projection_decompositions[0][1]));
6397 static_cast<double>(
6398 projection_decompositions[1][1]));
6401 static_cast<double>(
6402 projection_decompositions[2][1]));
6405 static_cast<double>(
6406 projection_decompositions[3][1]));
6410 static_cast<double>(
6411 projection_decompositions[0][1]));
6414 static_cast<double>(
6415 projection_decompositions[1][1]));
6418 static_cast<double>(
6419 projection_decompositions[2][1]));
6422 static_cast<double>(
6423 projection_decompositions[3][1]));
6428 x_dimension_perspective = x_max_perspective - x_min_perspective;
6429 y_dimension_perspective = y_max_perspective - y_min_perspective;
6431 std::multiset<SvgCell> cells;
6434 for (
const auto &patch : patches)
6436 n_subdivisions = patch.n_subdivisions;
6438 for (
unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
6440 for (
unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
6443 {compute_node(patch, i1, i2, 0, n_subdivisions),
6444 compute_node(patch, i1 + 1, i2, 0, n_subdivisions),
6445 compute_node(patch, i1, i2 + 1, 0, n_subdivisions),
6446 compute_node(patch, i1 + 1, i2 + 1, 0, n_subdivisions)}};
6449 patch.data.n_rows() == 0,
6452 patch.data.n_rows()));
6456 cell.vertices[0][0] = projected_vertices[0][0];
6457 cell.vertices[0][1] = projected_vertices[0][1];
6458 cell.vertices[0][2] = patch.data.n_rows() != 0 ?
6459 patch.data(0, i1 * d1 + i2 * d2) :
6462 cell.vertices[1][0] = projected_vertices[1][0];
6463 cell.vertices[1][1] = projected_vertices[1][1];
6464 cell.vertices[1][2] = patch.data.n_rows() != 0 ?
6465 patch.data(0, (i1 + 1) * d1 + i2 * d2) :
6468 cell.vertices[2][0] = projected_vertices[2][0];
6469 cell.vertices[2][1] = projected_vertices[2][1];
6470 cell.vertices[2][2] = patch.data.n_rows() != 0 ?
6471 patch.data(0, i1 * d1 + (i2 + 1) * d2) :
6474 cell.vertices[3][0] = projected_vertices[3][0];
6475 cell.vertices[3][1] = projected_vertices[3][1];
6476 cell.vertices[3][2] =
6477 patch.data.n_rows() != 0 ?
6478 patch.data(0, (i1 + 1) * d1 + (i2 + 1) * d2) :
6481 cell.projected_vertices[0] =
6482 svg_project_point(cell.vertices[0],
6487 cell.projected_vertices[1] =
6488 svg_project_point(cell.vertices[1],
6493 cell.projected_vertices[2] =
6494 svg_project_point(cell.vertices[2],
6499 cell.projected_vertices[3] =
6500 svg_project_point(cell.vertices[3],
6506 cell.center = .25 * (cell.vertices[0] + cell.vertices[1] +
6507 cell.vertices[2] + cell.vertices[3]);
6508 cell.projected_center = svg_project_point(cell.center,
6514 cell.depth = cell.center.distance(camera_position);
6524 width =
static_cast<unsigned int>(
6525 .5 + height * (x_dimension_perspective / y_dimension_perspective));
6526 unsigned int additional_width = 0;
6529 additional_width =
static_cast<unsigned int>(
6533 out <<
"<svg width=\"" << width + additional_width <<
"\" height=\"" 6534 << height <<
"\" xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\">" 6536 <<
" <rect width=\"" << width + additional_width <<
"\" height=\"" 6537 << height <<
"\" style=\"fill:white\"/>" <<
'\n' 6540 unsigned int triangle_counter = 0;
6543 for (
const auto &cell : cells)
6547 for (
unsigned int triangle_index = 0; triangle_index < 4;
6550 switch (triangle_index)
6553 points3d_triangle[0] = cell.vertices[0],
6554 points3d_triangle[1] = cell.vertices[1],
6555 points3d_triangle[2] = cell.center;
6558 points3d_triangle[0] = cell.vertices[1],
6559 points3d_triangle[1] = cell.vertices[3],
6560 points3d_triangle[2] = cell.center;
6563 points3d_triangle[0] = cell.vertices[3],
6564 points3d_triangle[1] = cell.vertices[2],
6565 points3d_triangle[2] = cell.center;
6568 points3d_triangle[0] = cell.vertices[2],
6569 points3d_triangle[1] = cell.vertices[0],
6570 points3d_triangle[2] = cell.center;
6577 svg_get_gradient_parameters(points3d_triangle);
6580 .667 - ((gradient_param[4] - z_min) / z_dimension) * .667;
6582 .667 - ((gradient_param[5] - z_min) / z_dimension) * .667;
6584 unsigned int start_r = 0;
6585 unsigned int start_g = 0;
6586 unsigned int start_b = 0;
6588 unsigned int stop_r = 0;
6589 unsigned int stop_g = 0;
6590 unsigned int stop_b = 0;
6592 unsigned int start_i =
static_cast<unsigned int>(start_h * 6.);
6593 unsigned int stop_i =
static_cast<unsigned int>(stop_h * 6.);
6595 double start_f = start_h * 6. - start_i;
6596 double start_q = 1. - start_f;
6598 double stop_f = stop_h * 6. - stop_i;
6599 double stop_q = 1. - stop_f;
6601 switch (start_i % 6)
6605 start_g =
static_cast<unsigned int>(.5 + 255. * start_f);
6608 start_r =
static_cast<unsigned int>(.5 + 255. * start_q),
6613 start_b =
static_cast<unsigned int>(.5 + 255. * start_f);
6616 start_g =
static_cast<unsigned int>(.5 + 255. * start_q),
6620 start_r =
static_cast<unsigned int>(.5 + 255. * start_f),
6625 start_b =
static_cast<unsigned int>(.5 + 255. * start_q);
6635 stop_g =
static_cast<unsigned int>(.5 + 255. * stop_f);
6638 stop_r =
static_cast<unsigned int>(.5 + 255. * stop_q),
6643 stop_b =
static_cast<unsigned int>(.5 + 255. * stop_f);
6646 stop_g =
static_cast<unsigned int>(.5 + 255. * stop_q),
6650 stop_r =
static_cast<unsigned int>(.5 + 255. * stop_f),
6655 stop_b =
static_cast<unsigned int>(.5 + 255. * stop_q);
6661 Point<3> gradient_start_point_3d, gradient_stop_point_3d;
6663 gradient_start_point_3d[0] = gradient_param[0];
6664 gradient_start_point_3d[1] = gradient_param[1];
6665 gradient_start_point_3d[2] = gradient_param[4];
6667 gradient_stop_point_3d[0] = gradient_param[2];
6668 gradient_stop_point_3d[1] = gradient_param[3];
6669 gradient_stop_point_3d[2] = gradient_param[5];
6672 svg_project_point(gradient_start_point_3d,
6678 svg_project_point(gradient_stop_point_3d,
6685 out <<
" <linearGradient id=\"" << triangle_counter
6686 <<
"\" gradientUnits=\"userSpaceOnUse\" " 6688 <<
static_cast<unsigned int>(
6690 ((gradient_start_point[0] - x_min_perspective) /
6691 x_dimension_perspective) *
6692 (width - (width / 100.) * 2. * margin_in_percent) +
6693 ((width / 100.) * margin_in_percent))
6696 <<
static_cast<unsigned int>(
6697 .5 + height - (height / 100.) * margin_in_percent -
6698 ((gradient_start_point[1] - y_min_perspective) /
6699 y_dimension_perspective) *
6700 (height - (height / 100.) * 2. * margin_in_percent))
6703 <<
static_cast<unsigned int>(
6705 ((gradient_stop_point[0] - x_min_perspective) /
6706 x_dimension_perspective) *
6707 (width - (width / 100.) * 2. * margin_in_percent) +
6708 ((width / 100.) * margin_in_percent))