Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
data_out_base.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1999 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16
19#include <deal.II/base/mpi.h>
24
26
27#include <algorithm>
28#include <cmath>
29#include <cstdint>
30#include <cstring>
31#include <ctime>
32#include <fstream>
33#include <iomanip>
34#include <limits>
35#include <memory>
36#include <set>
37#include <sstream>
38#include <vector>
39
40#ifdef DEAL_II_WITH_ZLIB
41# include <zlib.h>
42#endif
43
44#ifdef DEAL_II_WITH_HDF5
45# include <hdf5.h>
46#endif
47
48#include <boost/iostreams/copy.hpp>
49#include <boost/iostreams/device/back_inserter.hpp>
50#include <boost/iostreams/filtering_stream.hpp>
51#ifdef DEAL_II_WITH_ZLIB
52# include <boost/iostreams/filter/zlib.hpp>
53#endif
54
55
56
58
59// we need the following exception from a global function, so can't declare it
60// in the usual way inside a class
61namespace
62{
63 DeclException2(ExcUnexpectedInput,
64 std::string,
65 std::string,
66 << "Unexpected input: expected line\n <" << arg1
67 << ">\nbut got\n <" << arg2 << ">");
68}
69
70
71namespace
72{
73#ifdef DEAL_II_WITH_ZLIB
74 constexpr bool deal_ii_with_zlib = true;
75#else
76 constexpr bool deal_ii_with_zlib = false;
77#endif
78
79
80#ifdef DEAL_II_WITH_ZLIB
85 int
86 get_zlib_compression_level(const DataOutBase::CompressionLevel level)
87 {
88 switch (level)
89 {
91 return Z_NO_COMPRESSION;
93 return Z_BEST_SPEED;
95 return Z_BEST_COMPRESSION;
97 return Z_DEFAULT_COMPRESSION;
98 default:
99 Assert(false, ExcNotImplemented());
100 return Z_NO_COMPRESSION;
101 }
102 }
103
104# ifdef DEAL_II_WITH_MPI
109 int
110 get_boost_zlib_compression_level(const DataOutBase::CompressionLevel level)
111 {
112 switch (level)
113 {
115 return boost::iostreams::zlib::no_compression;
117 return boost::iostreams::zlib::best_speed;
119 return boost::iostreams::zlib::best_compression;
121 return boost::iostreams::zlib::default_compression;
122 default:
123 Assert(false, ExcNotImplemented());
124 return boost::iostreams::zlib::no_compression;
125 }
126 }
127# endif
128#endif
129
134 template <typename T>
135 std::string
136 compress_array(const std::vector<T> & data,
137 const DataOutBase::CompressionLevel compression_level)
138 {
139#ifdef DEAL_II_WITH_ZLIB
140 if (data.size() != 0)
141 {
142 const std::size_t uncompressed_size = (data.size() * sizeof(T));
143
144 // While zlib's compress2 uses unsigned long (which is 64bits
145 // on Linux), the vtu compression header stores the block size
146 // as an std::uint32_t (see below). While we could implement
147 // writing several smaller blocks, we haven't done that. Let's
148 // trigger an error for the user instead:
149 AssertThrow(uncompressed_size <=
150 std::numeric_limits<std::uint32_t>::max(),
152
153 // allocate a buffer for compressing data and do so
154 auto compressed_data_length = compressBound(uncompressed_size);
155 AssertThrow(compressed_data_length <=
156 std::numeric_limits<std::uint32_t>::max(),
158
159 std::vector<unsigned char> compressed_data(compressed_data_length);
160
161 int err = compress2(&compressed_data[0],
162 &compressed_data_length,
163 reinterpret_cast<const Bytef *>(data.data()),
164 uncompressed_size,
165 get_zlib_compression_level(compression_level));
166 (void)err;
167 Assert(err == Z_OK, ExcInternalError());
168
169 // Discard the unnecessary bytes
170 compressed_data.resize(compressed_data_length);
171
172 // now encode the compression header
173 const std::uint32_t compression_header[4] = {
174 1, /* number of blocks */
175 static_cast<std::uint32_t>(uncompressed_size), /* size of block */
176 static_cast<std::uint32_t>(
177 uncompressed_size), /* size of last block */
178 static_cast<std::uint32_t>(
179 compressed_data_length)}; /* list of compressed sizes of blocks */
180
181 const auto header_start =
182 reinterpret_cast<const unsigned char *>(&compression_header[0]);
183
185 {header_start, header_start + 4 * sizeof(std::uint32_t)}) +
186 Utilities::encode_base64(compressed_data));
187 }
188 else
189 return {};
190#else
191 (void)data;
192 (void)compression_level;
193 Assert(false,
194 ExcMessage("This function can only be called if cmake found "
195 "a working libz installation."));
196 return {};
197#endif
198 }
199
200
201
210 template <typename T>
211 std::string
212 vtu_stringize_array(const std::vector<T> & data,
213 const DataOutBase::CompressionLevel compression_level,
214 const int precision)
215 {
216 if (deal_ii_with_zlib &&
217 (compression_level != DataOutBase::CompressionLevel::plain_text))
218 {
219 // compress the data we have in memory
220 return compress_array(data, compression_level);
221 }
222 else
223 {
224 std::ostringstream stream;
225 stream.precision(precision);
226 for (const T &el : data)
227 stream << el << ' ';
228 return stream.str();
229 }
230 }
231
232
241 struct ParallelIntermediateHeader
242 {
243 std::uint64_t magic;
244 std::uint64_t version;
245 std::uint64_t compression;
246 std::uint64_t dimension;
247 std::uint64_t space_dimension;
248 std::uint64_t n_ranks;
249 std::uint64_t n_patches;
250 };
251} // namespace
252
253
254// some declarations of functions and locally used classes
255namespace DataOutBase
256{
257 namespace
258 {
264 class SvgCell
265 {
266 public:
267 // Center of the cell (three-dimensional)
269
274
279 float depth;
280
285
286 // Center of the cell (projected, two-dimensional)
288
292 bool
293 operator<(const SvgCell &) const;
294 };
295
296 bool
297 SvgCell::operator<(const SvgCell &e) const
298 {
299 // note the "wrong" order in which we sort the elements
300 return depth > e.depth;
301 }
302
303
304
310 class EpsCell2d
311 {
312 public:
317
323
328 float depth;
329
333 bool
334 operator<(const EpsCell2d &) const;
335 };
336
337 bool
338 EpsCell2d::operator<(const EpsCell2d &e) const
339 {
340 // note the "wrong" order in which we sort the elements
341 return depth > e.depth;
342 }
343
344
345
357 template <int dim, int spacedim, typename Number = double>
358 std::unique_ptr<Table<2, Number>>
359 create_global_data_table(const std::vector<Patch<dim, spacedim>> &patches)
361 // If there is nothing to write, just return
362 if (patches.size() == 0)
363 return std::make_unique<Table<2, Number>>();
364
365 // unlike in the main function, we don't have here the data_names field,
366 // so we initialize it with the number of data sets in the first patch.
367 // the equivalence of these two definitions is checked in the main
368 // function.
370 // we have to take care, however, whether the points are appended to the
371 // end of the patch.data table
372 const unsigned int n_data_sets = patches[0].points_are_available ?
373 (patches[0].data.n_rows() - spacedim) :
374 patches[0].data.n_rows();
375 const unsigned int n_data_points =
376 std::accumulate(patches.begin(),
377 patches.end(),
378 0U,
379 [](const unsigned int count,
380 const Patch<dim, spacedim> &patch) {
381 return count + patch.data.n_cols();
382 });
383
384 std::unique_ptr<Table<2, Number>> global_data_table =
385 std::make_unique<Table<2, Number>>(n_data_sets, n_data_points);
386
387 // loop over all patches
388 unsigned int next_value = 0;
389 for (const auto &patch : patches)
390 {
391 const unsigned int n_subdivisions = patch.n_subdivisions;
392 (void)n_subdivisions;
393
394 Assert((patch.data.n_rows() == n_data_sets &&
395 !patch.points_are_available) ||
396 (patch.data.n_rows() == n_data_sets + spacedim &&
397 patch.points_are_available),
398 ExcDimensionMismatch(patch.points_are_available ?
399 (n_data_sets + spacedim) :
400 n_data_sets,
401 patch.data.n_rows()));
402 Assert(patch.reference_cell != ReferenceCells::get_hypercube<dim>() ||
403 (n_data_sets == 0) ||
404 (patch.data.n_cols() ==
405 Utilities::fixed_power<dim>(n_subdivisions + 1)),
406 ExcInvalidDatasetSize(patch.data.n_cols(),
407 n_subdivisions + 1));
408
409 for (unsigned int i = 0; i < patch.data.n_cols(); ++i, ++next_value)
410 for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
411 (*global_data_table)[data_set][next_value] =
412 patch.data(data_set, i);
413 }
414 Assert(next_value == n_data_points, ExcInternalError());
415
416 return global_data_table;
417 }
418 } // namespace
419
420
421
423 : flags(false, true)
424 , node_dim(numbers::invalid_unsigned_int)
425 , num_cells(0)
426 {}
427
428
429
431 : flags(flags)
432 , node_dim(numbers::invalid_unsigned_int)
433 , num_cells(0)
434 {}
435
436
437
438 template <int dim>
439 void
440 DataOutFilter::write_point(const unsigned int index, const Point<dim> &p)
441 {
442 node_dim = dim;
443
444 Point<3> int_pt;
445 for (unsigned int d = 0; d < dim; ++d)
446 int_pt(d) = p(d);
447
448 const Map3DPoint::const_iterator it = existing_points.find(int_pt);
449 unsigned int internal_ind;
450
451 // If the point isn't in the set, or we're not filtering duplicate points,
452 // add it
454 {
455 internal_ind = existing_points.size();
456 existing_points.insert(std::make_pair(int_pt, internal_ind));
457 }
458 else
459 {
460 internal_ind = it->second;
461 }
462 // Now add the index to the list of filtered points
463 filtered_points[index] = internal_ind;
464 }
465
466
467
468 void
470 const unsigned int pt_index)
471 {
473
474 // (Re)-initialize counter at any first call to this method.
475 if (cell_index == 0)
476 num_cells = 1;
477 }
478
479
480
481 void
482 DataOutFilter::fill_node_data(std::vector<double> &node_data) const
483 {
484 node_data.resize(existing_points.size() * node_dim);
485
486 for (const auto &existing_point : existing_points)
487 {
488 for (unsigned int d = 0; d < node_dim; ++d)
489 node_data[node_dim * existing_point.second + d] =
490 existing_point.first(d);
491 }
492 }
493
494
495
496 void
497 DataOutFilter::fill_cell_data(const unsigned int local_node_offset,
498 std::vector<unsigned int> &cell_data) const
499 {
500 cell_data.resize(filtered_cells.size());
501
502 for (const auto &filtered_cell : filtered_cells)
503 {
504 cell_data[filtered_cell.first] =
505 filtered_cell.second + local_node_offset;
506 }
507 }
508
509
510
511 std::string
512 DataOutFilter::get_data_set_name(const unsigned int set_num) const
513 {
514 return data_set_names.at(set_num);
515 }
516
517
518
519 unsigned int
520 DataOutFilter::get_data_set_dim(const unsigned int set_num) const
521 {
522 return data_set_dims.at(set_num);
523 }
524
525
526
527 const double *
528 DataOutFilter::get_data_set(const unsigned int set_num) const
529 {
530 return data_sets[set_num].data();
531 }
532
533
534
535 unsigned int
537 {
538 return existing_points.size();
539 }
540
541
542
543 unsigned int
545 {
546 return num_cells;
547 }
548
549
550
551 unsigned int
553 {
554 return data_set_names.size();
555 }
556
557
558
559 void
561 {}
562
563
564
565 void
567 {}
568
569
570
571 template <int dim>
572 void
573 DataOutFilter::write_cell(const unsigned int index,
574 const unsigned int start,
575 const std::array<unsigned int, dim> &offsets)
576 {
577 ++num_cells;
578
579 const unsigned int base_entry =
581
582 switch (dim)
583 {
584 case 0:
585 {
586 internal_add_cell(base_entry + 0, start);
587 break;
588 }
589
590 case 1:
591 {
592 const unsigned int d1 = offsets[0];
593
594 internal_add_cell(base_entry + 0, start);
595 internal_add_cell(base_entry + 1, start + d1);
596 break;
597 }
598
599 case 2:
600 {
601 const unsigned int d1 = offsets[0];
602 const unsigned int d2 = offsets[1];
603
604 internal_add_cell(base_entry + 0, start);
605 internal_add_cell(base_entry + 1, start + d1);
606 internal_add_cell(base_entry + 2, start + d2 + d1);
607 internal_add_cell(base_entry + 3, start + d2);
608 break;
609 }
610
611 case 3:
612 {
613 const unsigned int d1 = offsets[0];
614 const unsigned int d2 = offsets[1];
615 const unsigned int d3 = offsets[2];
616
617 internal_add_cell(base_entry + 0, start);
618 internal_add_cell(base_entry + 1, start + d1);
619 internal_add_cell(base_entry + 2, start + d2 + d1);
620 internal_add_cell(base_entry + 3, start + d2);
621 internal_add_cell(base_entry + 4, start + d3);
622 internal_add_cell(base_entry + 5, start + d3 + d1);
623 internal_add_cell(base_entry + 6, start + d3 + d2 + d1);
624 internal_add_cell(base_entry + 7, start + d3 + d2);
625 break;
626 }
627
628 default:
629 Assert(false, ExcNotImplemented());
630 }
631 }
632
633
634
635 void
636 DataOutFilter::write_cell_single(const unsigned int index,
637 const unsigned int start,
638 const unsigned int n_points,
639 const ReferenceCell &reference_cell)
640 {
641 ++num_cells;
642
643 const unsigned int base_entry = index * n_points;
644
645 static const std::array<unsigned int, 5> table = {{0, 1, 3, 2, 4}};
646
647 for (unsigned int i = 0; i < n_points; ++i)
648 internal_add_cell(base_entry + i,
649 start + (reference_cell == ReferenceCells::Pyramid ?
650 table[i] :
651 i));
652 }
653
654
655
656 void
657 DataOutFilter::write_data_set(const std::string & name,
658 const unsigned int dimension,
659 const unsigned int set_num,
660 const Table<2, double> &data_vectors)
661 {
662 unsigned int new_dim;
663
664 // HDF5/XDMF output only supports 1d or 3d output, so force rearrangement if
665 // needed
666 if (flags.xdmf_hdf5_output && dimension != 1)
667 new_dim = 3;
668 else
669 new_dim = dimension;
670
671 // Record the data set name, dimension, and allocate space for it
672 data_set_names.push_back(name);
673 data_set_dims.push_back(new_dim);
674 data_sets.emplace_back(new_dim * existing_points.size());
675
676 // TODO: averaging, min/max, etc for merged vertices
677 for (unsigned int i = 0; i < filtered_points.size(); ++i)
678 {
679 const unsigned int r = filtered_points[i];
680
681 for (unsigned int d = 0; d < new_dim; ++d)
682 {
683 if (d < dimension)
684 data_sets.back()[r * new_dim + d] = data_vectors(set_num + d, i);
685 else
686 data_sets.back()[r * new_dim + d] = 0;
687 }
688 }
689 }
690} // namespace DataOutBase
691
692
693
694//----------------------------------------------------------------------//
695// Auxiliary data
696//----------------------------------------------------------------------//
697
698namespace
699{
700 const char *gmv_cell_type[4] = {"", "line 2", "quad 4", "hex 8"};
701
702 const char *ucd_cell_type[4] = {"pt", "line", "quad", "hex"};
703
704 const char *tecplot_cell_type[4] = {"", "lineseg", "quadrilateral", "brick"};
705
719 template <int dim, int spacedim>
720 std::array<unsigned int, 3>
721 extract_vtk_patch_info(const DataOutBase::Patch<dim, spacedim> &patch,
722 const bool write_higher_order_cells)
723 {
724 std::array<unsigned int, 3> vtk_cell_id = {
725 {/* cell type, tbd: */ numbers::invalid_unsigned_int,
726 /* # of cells, default: just one cell */ 1,
727 /* # of nodes, default: as many nodes as vertices */
728 patch.reference_cell.n_vertices()}};
729
730 if (write_higher_order_cells)
731 {
732 vtk_cell_id[0] = patch.reference_cell.vtk_lagrange_type();
733 vtk_cell_id[2] = patch.data.n_cols();
734 }
735 else if (patch.data.n_cols() == patch.reference_cell.n_vertices())
736 // One data set per vertex -> a linear cell
737 vtk_cell_id[0] = patch.reference_cell.vtk_linear_type();
738 else if (patch.reference_cell == ReferenceCells::Triangle &&
739 patch.data.n_cols() == 6)
740 {
742 vtk_cell_id[0] = patch.reference_cell.vtk_quadratic_type();
743 vtk_cell_id[2] = patch.data.n_cols();
744 }
746 patch.data.n_cols() == 10)
747 {
749 vtk_cell_id[0] = patch.reference_cell.vtk_quadratic_type();
750 vtk_cell_id[2] = patch.data.n_cols();
751 }
752 else if (patch.reference_cell.is_hyper_cube())
753 {
754 // For hypercubes, we support sub-divided linear cells
755 vtk_cell_id[0] = patch.reference_cell.vtk_linear_type();
756 vtk_cell_id[1] = Utilities::pow(patch.n_subdivisions, dim);
757 }
758 else
759 {
760 Assert(false, ExcNotImplemented());
761 }
762
763 return vtk_cell_id;
764 }
765
766 //----------------------------------------------------------------------//
767 // Auxiliary functions
768 //----------------------------------------------------------------------//
769
770 // For a given patch that corresponds to a hypercube cell, compute the
771 // location of a node interpolating the corner nodes linearly
772 // at the point lattice_location/n_subdivisions where lattice_location
773 // is a dim-dimensional integer vector. If the points are
774 // saved in the patch.data member, return the saved point instead.
775 template <int dim, int spacedim>
776 inline Point<spacedim>
777 get_equispaced_location(
779 const std::initializer_list<unsigned int> &lattice_location,
780 const unsigned int n_subdivisions)
781 {
782 // This function only makes sense when called on hypercube cells
784
785 Assert(lattice_location.size() == dim, ExcInternalError());
786
787 const unsigned int xstep = (dim > 0 ? *(lattice_location.begin() + 0) : 0);
788 const unsigned int ystep = (dim > 1 ? *(lattice_location.begin() + 1) : 0);
789 const unsigned int zstep = (dim > 2 ? *(lattice_location.begin() + 2) : 0);
790
791 // If the patch stores the locations of nodes (rather than of only the
792 // vertices), then obtain the location by direct lookup.
793 if (patch.points_are_available)
794 {
795 Assert(n_subdivisions == patch.n_subdivisions, ExcNotImplemented());
796
797 unsigned int point_no = 0;
798 switch (dim)
799 {
800 case 3:
801 AssertIndexRange(zstep, n_subdivisions + 1);
802 point_no += (n_subdivisions + 1) * (n_subdivisions + 1) * zstep;
804 case 2:
805 AssertIndexRange(ystep, n_subdivisions + 1);
806 point_no += (n_subdivisions + 1) * ystep;
808 case 1:
809 AssertIndexRange(xstep, n_subdivisions + 1);
810 point_no += xstep;
812 case 0:
813 // break here for dim<=3
814 break;
815
816 default:
817 Assert(false, ExcNotImplemented());
818 }
819 Point<spacedim> node;
820 for (unsigned int d = 0; d < spacedim; ++d)
821 node[d] = patch.data(patch.data.size(0) - spacedim + d, point_no);
822 return node;
823 }
824 else
825 // The patch does not store node locations, so we have to interpolate
826 // between its vertices:
827 {
828 if (dim == 0)
829 return patch.vertices[0];
830 else
831 {
832 // perform a dim-linear interpolation
833 const double stepsize = 1. / n_subdivisions;
834 const double xfrac = xstep * stepsize;
835
836 Point<spacedim> node =
837 (patch.vertices[1] * xfrac) + (patch.vertices[0] * (1 - xfrac));
838 if (dim > 1)
839 {
840 const double yfrac = ystep * stepsize;
841 node *= 1 - yfrac;
842 node += ((patch.vertices[3] * xfrac) +
843 (patch.vertices[2] * (1 - xfrac))) *
844 yfrac;
845 if (dim > 2)
846 {
847 const double zfrac = zstep * stepsize;
848 node *= (1 - zfrac);
849 node += (((patch.vertices[5] * xfrac) +
850 (patch.vertices[4] * (1 - xfrac))) *
851 (1 - yfrac) +
852 ((patch.vertices[7] * xfrac) +
853 (patch.vertices[6] * (1 - xfrac))) *
854 yfrac) *
855 zfrac;
856 }
857 }
858 return node;
859 }
860 }
861 }
862
863 // For a given patch, compute the nodes for arbitrary (non-hypercube) cells.
864 // If the points are saved in the patch.data member, return the saved point
865 // instead.
866 template <int dim, int spacedim>
867 inline Point<spacedim>
868 get_node_location(const DataOutBase::Patch<dim, spacedim> &patch,
869 const unsigned int node_index)
870 {
871 // Due to a historical accident, we are using a different indexing
872 // for pyramids in this file than we do where we create patches.
873 // So translate if necessary.
874 unsigned int point_no_actual = node_index;
876 {
878
879 static const std::array<unsigned int, 5> table = {{0, 1, 3, 2, 4}};
880 point_no_actual = table[node_index];
881 }
882
883 // If the patch stores the locations of nodes (rather than of only the
884 // vertices), then obtain the location by direct lookup.
885 if (patch.points_are_available)
886 {
887 Point<spacedim> node;
888 for (unsigned int d = 0; d < spacedim; ++d)
889 node[d] =
890 patch.data(patch.data.size(0) - spacedim + d, point_no_actual);
891 return node;
892 }
893 else
894 // The patch does not store node locations, so we have to interpolate
895 // between its vertices. This isn't currently implemented for anything
896 // other than one subdivision, but would go here.
897 //
898 // For n_subdivisions==1, the locations are simply those of vertices, so
899 // get the information from there.
900 {
902
903 return patch.vertices[point_no_actual];
904 }
905 }
906
907
908
914 template <int dim, int spacedim>
915 std::tuple<unsigned int, unsigned int>
916 count_nodes_and_cells(
917 const std::vector<DataOutBase::Patch<dim, spacedim>> &patches)
918 {
919 unsigned int n_nodes = 0;
920 unsigned int n_cells = 0;
921 for (const auto &patch : patches)
922 {
925 "The reference cell for this patch is set to 'Invalid', "
926 "but that is clearly not a valid choice. Did you forget "
927 "to set the reference cell for the patch?"));
928
929 if (patch.reference_cell.is_hyper_cube())
930 {
931 n_nodes += Utilities::fixed_power<dim>(patch.n_subdivisions + 1);
932 n_cells += Utilities::fixed_power<dim>(patch.n_subdivisions);
933 }
934 else
935 {
937 n_nodes += patch.reference_cell.n_vertices();
938 n_cells += 1;
939 }
940 }
941
942 return std::make_tuple(n_nodes, n_cells);
943 }
944
945
946
952 template <int dim, int spacedim>
953 std::tuple<unsigned int, unsigned int, unsigned int>
954 count_nodes_and_cells_and_points(
955 const std::vector<DataOutBase::Patch<dim, spacedim>> &patches,
956 const bool write_higher_order_cells)
957 {
958 unsigned int n_nodes = 0;
959 unsigned int n_cells = 0;
960 unsigned int n_points_and_n_cells = 0;
961
962 for (const auto &patch : patches)
963 {
964 if (patch.reference_cell.is_hyper_cube())
965 {
966 n_nodes += Utilities::fixed_power<dim>(patch.n_subdivisions + 1);
967
968 if (write_higher_order_cells)
969 {
970 // Write all of these nodes as a single higher-order cell. So
971 // add one to the number of cells, and update the number of
972 // points appropriately.
973 n_cells += 1;
974 n_points_and_n_cells +=
975 1 + Utilities::fixed_power<dim>(patch.n_subdivisions + 1);
976 }
977 else
978 {
979 // Write all of these nodes as a collection of d-linear
980 // cells. Add the number of sub-cells to the total number of
981 // cells, and then add one for each cell plus the number of
982 // vertices per cell for each subcell to the number of points.
983 const unsigned int n_subcells =
984 Utilities::fixed_power<dim>(patch.n_subdivisions);
985 n_cells += n_subcells;
986 n_points_and_n_cells +=
987 n_subcells * (1 + GeometryInfo<dim>::vertices_per_cell);
988 }
989 }
990 else
991 {
992 n_nodes += patch.data.n_cols();
993 n_cells += 1;
994 n_points_and_n_cells += patch.data.n_cols() + 1;
995 }
996 }
997
998 return std::make_tuple(n_nodes, n_cells, n_points_and_n_cells);
999 }
1000
1006 template <typename FlagsType>
1007 class StreamBase
1008 {
1009 public:
1010 /*
1011 * Constructor. Stores a reference to the output stream for immediate use.
1012 */
1013 StreamBase(std::ostream &stream, const FlagsType &flags)
1014 : selected_component(numbers::invalid_unsigned_int)
1015 , stream(stream)
1016 , flags(flags)
1017 {}
1018
1023 template <int dim>
1024 void
1025 write_point(const unsigned int, const Point<dim> &)
1026 {
1027 Assert(false,
1028 ExcMessage("The derived class you are using needs to "
1029 "reimplement this function if you want to call "
1030 "it."));
1031 }
1032
1038 void
1039 flush_points()
1040 {}
1041
1047 template <int dim>
1048 void
1049 write_cell(const unsigned int /*index*/,
1050 const unsigned int /*start*/,
1051 std::array<unsigned int, dim> & /*offsets*/)
1052 {
1053 Assert(false,
1054 ExcMessage("The derived class you are using needs to "
1055 "reimplement this function if you want to call "
1056 "it."));
1057 }
1058
1065 void
1066 write_cell_single(const unsigned int index,
1067 const unsigned int start,
1068 const unsigned int n_points,
1069 const ReferenceCell &reference_cell)
1070 {
1071 (void)index;
1072 (void)start;
1073 (void)n_points;
1074 (void)reference_cell;
1075
1076 Assert(false,
1077 ExcMessage("The derived class you are using needs to "
1078 "reimplement this function if you want to call "
1079 "it."));
1080 }
1081
1088 void
1089 flush_cells()
1090 {}
1091
1096 template <typename T>
1097 std::ostream &
1098 operator<<(const T &t)
1099 {
1100 stream << t;
1101 return stream;
1102 }
1103
1110 unsigned int selected_component;
1111
1112 protected:
1117 std::ostream &stream;
1118
1122 const FlagsType flags;
1123 };
1124
1128 class DXStream : public StreamBase<DataOutBase::DXFlags>
1129 {
1130 public:
1131 DXStream(std::ostream &stream, const DataOutBase::DXFlags &flags);
1132
1133 template <int dim>
1134 void
1135 write_point(const unsigned int index, const Point<dim> &);
1136
1145 template <int dim>
1146 void
1147 write_cell(const unsigned int index,
1148 const unsigned int start,
1149 const std::array<unsigned int, dim> &offsets);
1150
1157 template <typename data>
1158 void
1159 write_dataset(const unsigned int index, const std::vector<data> &values);
1160 };
1161
1165 class GmvStream : public StreamBase<DataOutBase::GmvFlags>
1166 {
1167 public:
1168 GmvStream(std::ostream &stream, const DataOutBase::GmvFlags &flags);
1169
1170 template <int dim>
1171 void
1172 write_point(const unsigned int index, const Point<dim> &);
1173
1182 template <int dim>
1183 void
1184 write_cell(const unsigned int index,
1185 const unsigned int start,
1186 const std::array<unsigned int, dim> &offsets);
1187 };
1188
1192 class TecplotStream : public StreamBase<DataOutBase::TecplotFlags>
1193 {
1194 public:
1195 TecplotStream(std::ostream &stream, const DataOutBase::TecplotFlags &flags);
1196
1197 template <int dim>
1198 void
1199 write_point(const unsigned int index, const Point<dim> &);
1200
1209 template <int dim>
1210 void
1211 write_cell(const unsigned int index,
1212 const unsigned int start,
1213 const std::array<unsigned int, dim> &offsets);
1214 };
1215
1219 class UcdStream : public StreamBase<DataOutBase::UcdFlags>
1220 {
1221 public:
1222 UcdStream(std::ostream &stream, const DataOutBase::UcdFlags &flags);
1223
1224 template <int dim>
1225 void
1226 write_point(const unsigned int index, const Point<dim> &);
1227
1238 template <int dim>
1239 void
1240 write_cell(const unsigned int index,
1241 const unsigned int start,
1242 const std::array<unsigned int, dim> &offsets);
1243
1250 template <typename data>
1251 void
1252 write_dataset(const unsigned int index, const std::vector<data> &values);
1253 };
1254
1258 class VtkStream : public StreamBase<DataOutBase::VtkFlags>
1259 {
1260 public:
1261 VtkStream(std::ostream &stream, const DataOutBase::VtkFlags &flags);
1262
1263 template <int dim>
1264 void
1265 write_point(const unsigned int index, const Point<dim> &);
1266
1275 template <int dim>
1276 void
1277 write_cell(const unsigned int index,
1278 const unsigned int start,
1279 const std::array<unsigned int, dim> &offsets);
1280
1284 void
1285 write_cell_single(const unsigned int index,
1286 const unsigned int start,
1287 const unsigned int n_points,
1288 const ReferenceCell &reference_cell);
1289
1297 template <int dim>
1298 void
1299 write_high_order_cell(const unsigned int start,
1300 const std::vector<unsigned> &connectivity);
1301 };
1302
1303
1304 //----------------------------------------------------------------------//
1305
1306 DXStream::DXStream(std::ostream &out, const DataOutBase::DXFlags &f)
1307 : StreamBase<DataOutBase::DXFlags>(out, f)
1308 {}
1309
1310
1311 template <int dim>
1312 void
1313 DXStream::write_point(const unsigned int, const Point<dim> &p)
1314 {
1315 if (flags.coordinates_binary)
1316 {
1317 float data[dim];
1318 for (unsigned int d = 0; d < dim; ++d)
1319 data[d] = p(d);
1320 stream.write(reinterpret_cast<const char *>(data), dim * sizeof(*data));
1321 }
1322 else
1323 {
1324 for (unsigned int d = 0; d < dim; ++d)
1325 stream << p(d) << '\t';
1326 stream << '\n';
1327 }
1328 }
1329
1330
1331
1332 // Separate these out to avoid an internal compiler error with intel 17
1334 {
1339 std::array<unsigned int, GeometryInfo<0>::vertices_per_cell>
1340 set_node_numbers(const unsigned int /*start*/,
1341 const std::array<unsigned int, 0> & /*d1*/)
1342 {
1343 Assert(false, ExcInternalError());
1344 return {};
1345 }
1346
1347
1348
1349 std::array<unsigned int, GeometryInfo<1>::vertices_per_cell>
1350 set_node_numbers(const unsigned int start,
1351 const std::array<unsigned int, 1> &offsets)
1352 {
1353 std::array<unsigned int, GeometryInfo<1>::vertices_per_cell> nodes;
1354 nodes[0] = start;
1355 nodes[1] = start + offsets[0];
1356 return nodes;
1357 }
1358
1359
1360
1361 std::array<unsigned int, GeometryInfo<2>::vertices_per_cell>
1362 set_node_numbers(const unsigned int start,
1363 const std::array<unsigned int, 2> &offsets)
1364
1365 {
1366 const unsigned int d1 = offsets[0];
1367 const unsigned int d2 = offsets[1];
1368
1369 std::array<unsigned int, GeometryInfo<2>::vertices_per_cell> nodes;
1370 nodes[0] = start;
1371 nodes[1] = start + d1;
1372 nodes[2] = start + d2;
1373 nodes[3] = start + d2 + d1;
1374 return nodes;
1375 }
1376
1377
1378
1379 std::array<unsigned int, GeometryInfo<3>::vertices_per_cell>
1380 set_node_numbers(const unsigned int start,
1381 const std::array<unsigned int, 3> &offsets)
1382 {
1383 const unsigned int d1 = offsets[0];
1384 const unsigned int d2 = offsets[1];
1385 const unsigned int d3 = offsets[2];
1386
1387 std::array<unsigned int, GeometryInfo<3>::vertices_per_cell> nodes;
1388 nodes[0] = start;
1389 nodes[1] = start + d1;
1390 nodes[2] = start + d2;
1391 nodes[3] = start + d2 + d1;
1392 nodes[4] = start + d3;
1393 nodes[5] = start + d3 + d1;
1394 nodes[6] = start + d3 + d2;
1395 nodes[7] = start + d3 + d2 + d1;
1396 return nodes;
1397 }
1398 } // namespace DataOutBaseImplementation
1399
1400
1401
1402 template <int dim>
1403 void
1404 DXStream::write_cell(const unsigned int,
1405 const unsigned int start,
1406 const std::array<unsigned int, dim> &offsets)
1407 {
1408 const auto nodes =
1409 DataOutBaseImplementation::set_node_numbers(start, offsets);
1410
1411 if (flags.int_binary)
1412 {
1413 std::array<unsigned int, GeometryInfo<dim>::vertices_per_cell> temp;
1414 for (unsigned int i = 0; i < nodes.size(); ++i)
1415 temp[i] = nodes[GeometryInfo<dim>::dx_to_deal[i]];
1416 stream.write(reinterpret_cast<const char *>(temp.data()),
1417 temp.size() * sizeof(temp[0]));
1418 }
1419 else
1420 {
1421 for (unsigned int i = 0; i < nodes.size() - 1; ++i)
1422 stream << nodes[GeometryInfo<dim>::dx_to_deal[i]] << '\t';
1423 stream << nodes[GeometryInfo<dim>::dx_to_deal[nodes.size() - 1]]
1424 << '\n';
1425 }
1426 }
1427
1428
1429
1430 template <typename data>
1431 inline void
1432 DXStream::write_dataset(const unsigned int, const std::vector<data> &values)
1433 {
1434 if (flags.data_binary)
1435 {
1436 stream.write(reinterpret_cast<const char *>(values.data()),
1437 values.size() * sizeof(data));
1438 }
1439 else
1440 {
1441 for (unsigned int i = 0; i < values.size(); ++i)
1442 stream << '\t' << values[i];
1443 stream << '\n';
1444 }
1445 }
1446
1447
1448
1449 //----------------------------------------------------------------------//
1450
1451 GmvStream::GmvStream(std::ostream &out, const DataOutBase::GmvFlags &f)
1452 : StreamBase<DataOutBase::GmvFlags>(out, f)
1453 {}
1454
1455
1456 template <int dim>
1457 void
1458 GmvStream::write_point(const unsigned int, const Point<dim> &p)
1459 {
1460 Assert(selected_component != numbers::invalid_unsigned_int,
1462 stream << p(selected_component) << ' ';
1463 }
1464
1465
1466
1467 template <int dim>
1468 void
1469 GmvStream::write_cell(const unsigned int,
1470 const unsigned int s,
1471 const std::array<unsigned int, dim> &offsets)
1472 {
1473 // Vertices are numbered starting with one.
1474 const unsigned int start = s + 1;
1475 stream << gmv_cell_type[dim] << '\n';
1476
1477 switch (dim)
1478 {
1479 case 0:
1480 {
1481 stream << start;
1482 break;
1483 }
1484
1485 case 1:
1486 {
1487 const unsigned int d1 = offsets[0];
1488 stream << start;
1489 stream << '\t' << start + d1;
1490 break;
1491 }
1492
1493 case 2:
1494 {
1495 const unsigned int d1 = offsets[0];
1496 const unsigned int d2 = offsets[1];
1497 stream << start;
1498 stream << '\t' << start + d1;
1499 stream << '\t' << start + d2 + d1 << '\t' << start + d2;
1500 break;
1501 }
1502
1503 case 3:
1504 {
1505 const unsigned int d1 = offsets[0];
1506 const unsigned int d2 = offsets[1];
1507 const unsigned int d3 = offsets[2];
1508 stream << start;
1509 stream << '\t' << start + d1;
1510 stream << '\t' << start + d2 + d1 << '\t' << start + d2;
1511 stream << '\t' << start + d3 << '\t' << start + d3 + d1 << '\t'
1512 << start + d3 + d2 + d1 << '\t' << start + d3 + d2;
1513 break;
1514 }
1515
1516 default:
1517 Assert(false, ExcNotImplemented());
1518 }
1519 stream << '\n';
1520 }
1521
1522
1523
1524 TecplotStream::TecplotStream(std::ostream & out,
1526 : StreamBase<DataOutBase::TecplotFlags>(out, f)
1527 {}
1528
1529
1530 template <int dim>
1531 void
1532 TecplotStream::write_point(const unsigned int, const Point<dim> &p)
1533 {
1534 Assert(selected_component != numbers::invalid_unsigned_int,
1536 stream << p(selected_component) << '\n';
1537 }
1538
1539
1540
1541 template <int dim>
1542 void
1543 TecplotStream::write_cell(const unsigned int,
1544 const unsigned int s,
1545 const std::array<unsigned int, dim> &offsets)
1546 {
1547 const unsigned int start = s + 1;
1548
1549 switch (dim)
1550 {
1551 case 0:
1552 {
1553 stream << start;
1554 break;
1555 }
1556
1557 case 1:
1558 {
1559 const unsigned int d1 = offsets[0];
1560 stream << start;
1561 stream << '\t' << start + d1;
1562 break;
1563 }
1564
1565 case 2:
1566 {
1567 const unsigned int d1 = offsets[0];
1568 const unsigned int d2 = offsets[1];
1569 stream << start;
1570 stream << '\t' << start + d1;
1571 stream << '\t' << start + d2 + d1 << '\t' << start + d2;
1572 break;
1573 }
1574
1575 case 3:
1576 {
1577 const unsigned int d1 = offsets[0];
1578 const unsigned int d2 = offsets[1];
1579 const unsigned int d3 = offsets[2];
1580 stream << start;
1581 stream << '\t' << start + d1;
1582 stream << '\t' << start + d2 + d1 << '\t' << start + d2;
1583 stream << '\t' << start + d3 << '\t' << start + d3 + d1 << '\t'
1584 << start + d3 + d2 + d1 << '\t' << start + d3 + d2;
1585 break;
1586 }
1587
1588 default:
1589 Assert(false, ExcNotImplemented());
1590 }
1591 stream << '\n';
1592 }
1593
1594
1595
1596 UcdStream::UcdStream(std::ostream &out, const DataOutBase::UcdFlags &f)
1597 : StreamBase<DataOutBase::UcdFlags>(out, f)
1598 {}
1599
1600
1601 template <int dim>
1602 void
1603 UcdStream::write_point(const unsigned int index, const Point<dim> &p)
1604 {
1605 stream << index + 1 << " ";
1606 // write out coordinates
1607 for (unsigned int i = 0; i < dim; ++i)
1608 stream << p(i) << ' ';
1609 // fill with zeroes
1610 for (unsigned int i = dim; i < 3; ++i)
1611 stream << "0 ";
1612 stream << '\n';
1613 }
1614
1615
1616
1617 template <int dim>
1618 void
1619 UcdStream::write_cell(const unsigned int index,
1620 const unsigned int start,
1621 const std::array<unsigned int, dim> &offsets)
1622 {
1623 const auto nodes =
1624 DataOutBaseImplementation::set_node_numbers(start, offsets);
1625
1626 // Write out all cells and remember that all indices must be shifted by one.
1627 stream << index + 1 << "\t0 " << ucd_cell_type[dim];
1628 for (unsigned int i = 0; i < nodes.size(); ++i)
1629 stream << '\t' << nodes[GeometryInfo<dim>::ucd_to_deal[i]] + 1;
1630 stream << '\n';
1631 }
1632
1633
1634
1635 template <typename data>
1636 inline void
1637 UcdStream::write_dataset(const unsigned int index,
1638 const std::vector<data> &values)
1639 {
1640 stream << index + 1;
1641 for (unsigned int i = 0; i < values.size(); ++i)
1642 stream << '\t' << values[i];
1643 stream << '\n';
1644 }
1645
1646
1647
1648 //----------------------------------------------------------------------//
1649
1650 VtkStream::VtkStream(std::ostream &out, const DataOutBase::VtkFlags &f)
1651 : StreamBase<DataOutBase::VtkFlags>(out, f)
1652 {}
1653
1654
1655 template <int dim>
1656 void
1657 VtkStream::write_point(const unsigned int, const Point<dim> &p)
1658 {
1659 // write out coordinates
1660 stream << p;
1661 // fill with zeroes
1662 for (unsigned int i = dim; i < 3; ++i)
1663 stream << " 0";
1664 stream << '\n';
1665 }
1666
1667
1668
1669 template <int dim>
1670 void
1671 VtkStream::write_cell(const unsigned int,
1672 const unsigned int start,
1673 const std::array<unsigned int, dim> &offsets)
1674 {
1675 stream << GeometryInfo<dim>::vertices_per_cell << '\t';
1676
1677 switch (dim)
1678 {
1679 case 0:
1680 {
1681 stream << start;
1682 break;
1683 }
1684
1685 case 1:
1686 {
1687 const unsigned int d1 = offsets[0];
1688 stream << start;
1689 stream << '\t' << start + d1;
1690 break;
1691 }
1692
1693 case 2:
1694 {
1695 const unsigned int d1 = offsets[0];
1696 const unsigned int d2 = offsets[1];
1697 stream << start;
1698 stream << '\t' << start + d1;
1699 stream << '\t' << start + d2 + d1 << '\t' << start + d2;
1700 break;
1701 }
1702
1703 case 3:
1704 {
1705 const unsigned int d1 = offsets[0];
1706 const unsigned int d2 = offsets[1];
1707 const unsigned int d3 = offsets[2];
1708 stream << start;
1709 stream << '\t' << start + d1;
1710 stream << '\t' << start + d2 + d1 << '\t' << start + d2;
1711 stream << '\t' << start + d3 << '\t' << start + d3 + d1 << '\t'
1712 << start + d3 + d2 + d1 << '\t' << start + d3 + d2;
1713 break;
1714 }
1715
1716 default:
1717 Assert(false, ExcNotImplemented());
1718 }
1719 stream << '\n';
1720 }
1721
1722
1723
1724 void
1725 VtkStream::write_cell_single(const unsigned int index,
1726 const unsigned int start,
1727 const unsigned int n_points,
1728 const ReferenceCell &reference_cell)
1729 {
1730 (void)index;
1731
1732 static const std::array<unsigned int, 5> table = {{0, 1, 3, 2, 4}};
1733
1734 stream << '\t' << n_points;
1735 for (unsigned int i = 0; i < n_points; ++i)
1736 stream << '\t'
1737 << start +
1738 (reference_cell == ReferenceCells::Pyramid ? table[i] : i);
1739 stream << '\n';
1740 }
1741
1742 template <int dim>
1743 void
1744 VtkStream::write_high_order_cell(const unsigned int start,
1745 const std::vector<unsigned> &connectivity)
1746 {
1747 stream << connectivity.size();
1748 for (const auto &c : connectivity)
1749 stream << '\t' << start + c;
1750 stream << '\n';
1751 }
1752} // namespace
1753
1754
1755
1756namespace DataOutBase
1757{
1758 const unsigned int Deal_II_IntermediateFlags::format_version = 4;
1759
1760
1761 template <int dim, int spacedim>
1762 const unsigned int Patch<dim, spacedim>::space_dim;
1763
1764
1765 template <int dim, int spacedim>
1766 const unsigned int Patch<dim, spacedim>::no_neighbor;
1767
1768
1769 template <int dim, int spacedim>
1771 : patch_index(no_neighbor)
1772 , n_subdivisions(1)
1773 , points_are_available(false)
1774 , reference_cell(ReferenceCells::Invalid)
1775 // all the other data has a constructor of its own, except for the "neighbors"
1776 // field, which we set to invalid values.
1777 {
1778 for (const unsigned int i : GeometryInfo<dim>::face_indices())
1780
1781 AssertIndexRange(dim, spacedim + 1);
1782 Assert(spacedim <= 3, ExcNotImplemented());
1783 }
1784
1785
1786
1787 template <int dim, int spacedim>
1788 bool
1790 {
1791 if (reference_cell != patch.reference_cell)
1792 return false;
1793
1794 // TODO: make tolerance relative
1795 const double epsilon = 3e-16;
1796 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
1797 if (vertices[i].distance(patch.vertices[i]) > epsilon)
1798 return false;
1799
1800 for (const unsigned int i : GeometryInfo<dim>::face_indices())
1801 if (neighbors[i] != patch.neighbors[i])
1802 return false;
1803
1804 if (patch_index != patch.patch_index)
1805 return false;
1806
1807 if (n_subdivisions != patch.n_subdivisions)
1808 return false;
1809
1810 if (points_are_available != patch.points_are_available)
1811 return false;
1812
1813 if (data.n_rows() != patch.data.n_rows())
1814 return false;
1815
1816 if (data.n_cols() != patch.data.n_cols())
1817 return false;
1818
1819 for (unsigned int i = 0; i < data.n_rows(); ++i)
1820 for (unsigned int j = 0; j < data.n_cols(); ++j)
1821 if (data[i][j] != patch.data[i][j])
1822 return false;
1823
1824 return true;
1825 }
1826
1827
1828
1829 template <int dim, int spacedim>
1830 std::size_t
1832 {
1833 return (sizeof(vertices) / sizeof(vertices[0]) *
1835 sizeof(neighbors) / sizeof(neighbors[0]) *
1840 MemoryConsumption::memory_consumption(points_are_available) +
1841 sizeof(reference_cell));
1842 }
1843
1844
1845
1846 template <int dim, int spacedim>
1847 void
1849 {
1850 std::swap(vertices, other_patch.vertices);
1851 std::swap(neighbors, other_patch.neighbors);
1852 std::swap(patch_index, other_patch.patch_index);
1853 std::swap(n_subdivisions, other_patch.n_subdivisions);
1854 data.swap(other_patch.data);
1855 std::swap(points_are_available, other_patch.points_are_available);
1856 std::swap(reference_cell, other_patch.reference_cell);
1857 }
1858
1859
1860
1861 template <int spacedim>
1862 const unsigned int Patch<0, spacedim>::space_dim;
1863
1864
1865 template <int spacedim>
1866 const unsigned int Patch<0, spacedim>::no_neighbor;
1867
1868
1869 template <int spacedim>
1870 unsigned int Patch<0, spacedim>::neighbors[1] = {
1872
1873 template <int spacedim>
1874 const unsigned int Patch<0, spacedim>::n_subdivisions = 1;
1875
1876 template <int spacedim>
1879
1880 template <int spacedim>
1882 : patch_index(no_neighbor)
1883 , points_are_available(false)
1884 {
1885 Assert(spacedim <= 3, ExcNotImplemented());
1886 }
1887
1888
1889
1890 template <int spacedim>
1891 bool
1893 {
1894 const unsigned int dim = 0;
1895
1896 // TODO: make tolerance relative
1897 const double epsilon = 3e-16;
1898 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
1899 if (vertices[i].distance(patch.vertices[i]) > epsilon)
1900 return false;
1901
1902 if (patch_index != patch.patch_index)
1903 return false;
1904
1905 if (points_are_available != patch.points_are_available)
1906 return false;
1907
1908 if (data.n_rows() != patch.data.n_rows())
1909 return false;
1910
1911 if (data.n_cols() != patch.data.n_cols())
1912 return false;
1913
1914 for (unsigned int i = 0; i < data.n_rows(); ++i)
1915 for (unsigned int j = 0; j < data.n_cols(); ++j)
1916 if (data[i][j] != patch.data[i][j])
1917 return false;
1918
1919 return true;
1920 }
1921
1922
1923
1924 template <int spacedim>
1925 std::size_t
1927 {
1928 return (sizeof(vertices) / sizeof(vertices[0]) *
1931 MemoryConsumption::memory_consumption(points_are_available));
1932 }
1933
1934
1935
1936 template <int spacedim>
1937 void
1939 {
1940 std::swap(vertices, other_patch.vertices);
1941 std::swap(patch_index, other_patch.patch_index);
1942 data.swap(other_patch.data);
1943 std::swap(points_are_available, other_patch.points_are_available);
1944 }
1945
1946
1947
1948 UcdFlags::UcdFlags(const bool write_preamble)
1949 : write_preamble(write_preamble)
1950 {}
1951
1952
1953
1955 {
1956 space_dimension_labels.emplace_back("x");
1957 space_dimension_labels.emplace_back("y");
1958 space_dimension_labels.emplace_back("z");
1959 }
1960
1961
1962
1963 GnuplotFlags::GnuplotFlags(const std::vector<std::string> &labels)
1964 : space_dimension_labels(labels)
1965 {}
1966
1967
1968
1969 std::size_t
1971 {
1973 }
1974
1975
1976
1977 PovrayFlags::PovrayFlags(const bool smooth,
1978 const bool bicubic_patch,
1979 const bool external_data)
1980 : smooth(smooth)
1981 , bicubic_patch(bicubic_patch)
1982 , external_data(external_data)
1983 {}
1984
1985
1986 DataOutFilterFlags::DataOutFilterFlags(const bool filter_duplicate_vertices,
1987 const bool xdmf_hdf5_output)
1988 : filter_duplicate_vertices(filter_duplicate_vertices)
1989 , xdmf_hdf5_output(xdmf_hdf5_output)
1990 {}
1991
1992
1993 void
1995 {
1996 prm.declare_entry(
1997 "Filter duplicate vertices",
1998 "false",
2000 "Whether to remove duplicate vertex values. deal.II duplicates "
2001 "vertices once for each adjacent cell so that it can output "
2002 "discontinuous quantities for which there may be more than one "
2003 "value for each vertex position. Setting this flag to "
2004 "'true' will merge all of these values by selecting a "
2005 "random one and outputting this as 'the' value for the vertex. "
2006 "As long as the data to be output corresponds to continuous "
2007 "fields, merging vertices has no effect. On the other hand, "
2008 "if the data to be output corresponds to discontinuous fields "
2009 "(either because you are using a discontinuous finite element, "
2010 "or because you are using a DataPostprocessor that yields "
2011 "discontinuous data, or because the data to be output has been "
2012 "produced by entirely different means), then the data in the "
2013 "output file no longer faithfully represents the underlying data "
2014 "because the discontinuous field has been replaced by a "
2015 "continuous one. Note also that the filtering can not occur "
2016 "on processor boundaries. Thus, a filtered discontinuous field "
2017 "looks like a continuous field inside of a subdomain, "
2018 "but like a discontinuous field at the subdomain boundary."
2019 "\n\n"
2020 "In any case, filtering results in drastically smaller output "
2021 "files (smaller by about a factor of 2^dim).");
2022 prm.declare_entry(
2023 "XDMF HDF5 output",
2024 "false",
2026 "Whether the data will be used in an XDMF/HDF5 combination.");
2027 }
2028
2029
2030
2031 void
2033 {
2034 filter_duplicate_vertices = prm.get_bool("Filter duplicate vertices");
2035 xdmf_hdf5_output = prm.get_bool("XDMF HDF5 output");
2036 }
2037
2038
2039
2040 DXFlags::DXFlags(const bool write_neighbors,
2041 const bool int_binary,
2042 const bool coordinates_binary,
2043 const bool data_binary)
2044 : write_neighbors(write_neighbors)
2045 , int_binary(int_binary)
2046 , coordinates_binary(coordinates_binary)
2047 , data_binary(data_binary)
2048 , data_double(false)
2049 {}
2050
2051
2052 void
2054 {
2055 prm.declare_entry("Write neighbors",
2056 "true",
2058 "A boolean field indicating whether neighborship "
2059 "information between cells is to be written to the "
2060 "OpenDX output file");
2061 prm.declare_entry("Integer format",
2062 "ascii",
2063 Patterns::Selection("ascii|32|64"),
2064 "Output format of integer numbers, which is "
2065 "either a text representation (ascii) or binary integer "
2066 "values of 32 or 64 bits length");
2067 prm.declare_entry("Coordinates format",
2068 "ascii",
2069 Patterns::Selection("ascii|32|64"),
2070 "Output format of vertex coordinates, which is "
2071 "either a text representation (ascii) or binary "
2072 "floating point values of 32 or 64 bits length");
2073 prm.declare_entry("Data format",
2074 "ascii",
2075 Patterns::Selection("ascii|32|64"),
2076 "Output format of data values, which is "
2077 "either a text representation (ascii) or binary "
2078 "floating point values of 32 or 64 bits length");
2079 }
2080
2081
2082
2083 void
2085 {
2086 write_neighbors = prm.get_bool("Write neighbors");
2087 // TODO:[GK] Read the new parameters
2088 }
2089
2090
2091
2092 void
2094 {
2095 prm.declare_entry("Write preamble",
2096 "true",
2098 "A flag indicating whether a comment should be "
2099 "written to the beginning of the output file "
2100 "indicating date and time of creation as well "
2101 "as the creating program");
2102 }
2103
2104
2105
2106 void
2108 {
2109 write_preamble = prm.get_bool("Write preamble");
2110 }
2111
2112
2113
2114 SvgFlags::SvgFlags(const unsigned int height_vector,
2115 const int azimuth_angle,
2116 const int polar_angle,
2117 const unsigned int line_thickness,
2118 const bool margin,
2119 const bool draw_colorbar)
2120 : height(4000)
2121 , width(0)
2122 , height_vector(height_vector)
2123 , azimuth_angle(azimuth_angle)
2124 , polar_angle(polar_angle)
2125 , line_thickness(line_thickness)
2126 , margin(margin)
2127 , draw_colorbar(draw_colorbar)
2128 {}
2129
2130
2131
2132 void
2134 {
2135 prm.declare_entry("Use smooth triangles",
2136 "false",
2138 "A flag indicating whether POVRAY should use smoothed "
2139 "triangles instead of the usual ones");
2140 prm.declare_entry("Use bicubic patches",
2141 "false",
2143 "Whether POVRAY should use bicubic patches");
2144 prm.declare_entry("Include external file",
2145 "true",
2147 "Whether camera and lighting information should "
2148 "be put into an external file \"data.inc\" or into "
2149 "the POVRAY input file");
2150 }
2151
2152
2153
2154 void
2156 {
2157 smooth = prm.get_bool("Use smooth triangles");
2158 bicubic_patch = prm.get_bool("Use bicubic patches");
2159 external_data = prm.get_bool("Include external file");
2160 }
2161
2162
2163
2164 EpsFlags::EpsFlags(const unsigned int height_vector,
2165 const unsigned int color_vector,
2166 const SizeType size_type,
2167 const unsigned int size,
2168 const double line_width,
2169 const double azimut_angle,
2170 const double turn_angle,
2171 const double z_scaling,
2172 const bool draw_mesh,
2173 const bool draw_cells,
2174 const bool shade_cells,
2175 const ColorFunction color_function)
2176 : height_vector(height_vector)
2177 , color_vector(color_vector)
2178 , size_type(size_type)
2179 , size(size)
2180 , line_width(line_width)
2181 , azimut_angle(azimut_angle)
2182 , turn_angle(turn_angle)
2183 , z_scaling(z_scaling)
2184 , draw_mesh(draw_mesh)
2185 , draw_cells(draw_cells)
2186 , shade_cells(shade_cells)
2187 , color_function(color_function)
2188 {}
2189
2190
2191
2194 const double xmin,
2195 const double xmax)
2196 {
2197 RgbValues rgb_values = {0, 0, 0};
2198
2199 // A difficult color scale:
2200 // xmin = black (1)
2201 // 3/4*xmin+1/4*xmax = blue (2)
2202 // 1/2*xmin+1/2*xmax = green (3)
2203 // 1/4*xmin+3/4*xmax = red (4)
2204 // xmax = white (5)
2205 // Makes the following color functions:
2206 //
2207 // red green blue
2208 // __
2209 // / /\ / /\ /
2210 // ____/ __/ \/ / \__/
2211
2212 // { 0 (1) - (3)
2213 // r = { ( 4*x-2*xmin+2*xmax)/(xmax-xmin) (3) - (4)
2214 // { 1 (4) - (5)
2215 //
2216 // { 0 (1) - (2)
2217 // g = { ( 4*x-3*xmin- xmax)/(xmax-xmin) (2) - (3)
2218 // { (-4*x+ xmin+3*xmax)/(xmax-xmin) (3) - (4)
2219 // { ( 4*x- xmin-3*xmax)/(xmax-xmin) (4) - (5)
2220 //
2221 // { ( 4*x-4*xmin )/(xmax-xmin) (1) - (2)
2222 // b = { (-4*x+2*xmin+2*xmax)/(xmax-xmin) (2) - (3)
2223 // { 0 (3) - (4)
2224 // { ( 4*x- xmin-3*xmax)/(xmax-xmin) (4) - (5)
2225
2226 double sum = xmax + xmin;
2227 double sum13 = xmin + 3 * xmax;
2228 double sum22 = 2 * xmin + 2 * xmax;
2229 double sum31 = 3 * xmin + xmax;
2230 double dif = xmax - xmin;
2231 double rezdif = 1.0 / dif;
2232
2233 int where;
2234
2235 if (x < (sum31) / 4)
2236 where = 0;
2237 else if (x < (sum22) / 4)
2238 where = 1;
2239 else if (x < (sum13) / 4)
2240 where = 2;
2241 else
2242 where = 3;
2243
2244 if (dif != 0)
2245 {
2246 switch (where)
2247 {
2248 case 0:
2249 rgb_values.red = 0;
2250 rgb_values.green = 0;
2251 rgb_values.blue = (x - xmin) * 4. * rezdif;
2252 break;
2253 case 1:
2254 rgb_values.red = 0;
2255 rgb_values.green = (4 * x - 3 * xmin - xmax) * rezdif;
2256 rgb_values.blue = (sum22 - 4. * x) * rezdif;
2257 break;
2258 case 2:
2259 rgb_values.red = (4 * x - 2 * sum) * rezdif;
2260 rgb_values.green = (xmin + 3 * xmax - 4 * x) * rezdif;
2261 rgb_values.blue = 0;
2262 break;
2263 case 3:
2264 rgb_values.red = 1;
2265 rgb_values.green = (4 * x - xmin - 3 * xmax) * rezdif;
2266 rgb_values.blue = (4. * x - sum13) * rezdif;
2267 break;
2268 default:
2269 break;
2270 }
2271 }
2272 else // White
2273 rgb_values.red = rgb_values.green = rgb_values.blue = 1;
2274
2275 return rgb_values;
2276 }
2277
2278
2279
2282 const double xmin,
2283 const double xmax)
2284 {
2285 EpsFlags::RgbValues rgb_values;
2286 rgb_values.red = rgb_values.blue = rgb_values.green =
2287 (x - xmin) / (xmax - xmin);
2288 return rgb_values;
2289 }
2290
2291
2292
2295 const double xmin,
2296 const double xmax)
2297 {
2298 EpsFlags::RgbValues rgb_values;
2299 rgb_values.red = rgb_values.blue = rgb_values.green =
2300 1 - (x - xmin) / (xmax - xmin);
2301 return rgb_values;
2302 }
2303
2304
2305
2306 void
2308 {
2309 prm.declare_entry("Index of vector for height",
2310 "0",
2312 "Number of the input vector that is to be used to "
2313 "generate height information");
2314 prm.declare_entry("Index of vector for color",
2315 "0",
2317 "Number of the input vector that is to be used to "
2318 "generate color information");
2319 prm.declare_entry("Scale to width or height",
2320 "width",
2321 Patterns::Selection("width|height"),
2322 "Whether width or height should be scaled to match "
2323 "the given size");
2324 prm.declare_entry("Size (width or height) in eps units",
2325 "300",
2327 "The size (width or height) to which the eps output "
2328 "file is to be scaled");
2329 prm.declare_entry("Line widths in eps units",
2330 "0.5",
2332 "The width in which the postscript renderer is to "
2333 "plot lines");
2334 prm.declare_entry("Azimut angle",
2335 "60",
2336 Patterns::Double(0, 180),
2337 "Angle of the viewing position against the vertical "
2338 "axis");
2339 prm.declare_entry("Turn angle",
2340 "30",
2341 Patterns::Double(0, 360),
2342 "Angle of the viewing direction against the y-axis");
2343 prm.declare_entry("Scaling for z-axis",
2344 "1",
2346 "Scaling for the z-direction relative to the scaling "
2347 "used in x- and y-directions");
2348 prm.declare_entry("Draw mesh lines",
2349 "true",
2351 "Whether the mesh lines, or only the surface should be "
2352 "drawn");
2353 prm.declare_entry("Fill interior of cells",
2354 "true",
2356 "Whether only the mesh lines, or also the interior of "
2357 "cells should be plotted. If this flag is false, then "
2358 "one can see through the mesh");
2359 prm.declare_entry("Color shading of interior of cells",
2360 "true",
2362 "Whether the interior of cells shall be shaded");
2363 prm.declare_entry("Color function",
2364 "default",
2366 "default|grey scale|reverse grey scale"),
2367 "Name of a color function used to colorize mesh lines "
2368 "and/or cell interiors");
2369 }
2370
2371
2372
2373 void
2375 {
2376 height_vector = prm.get_integer("Index of vector for height");
2377 color_vector = prm.get_integer("Index of vector for color");
2378 if (prm.get("Scale to width or height") == "width")
2379 size_type = width;
2380 else
2381 size_type = height;
2382 size = prm.get_integer("Size (width or height) in eps units");
2383 line_width = prm.get_double("Line widths in eps units");
2384 azimut_angle = prm.get_double("Azimut angle");
2385 turn_angle = prm.get_double("Turn angle");
2386 z_scaling = prm.get_double("Scaling for z-axis");
2387 draw_mesh = prm.get_bool("Draw mesh lines");
2388 draw_cells = prm.get_bool("Fill interior of cells");
2389 shade_cells = prm.get_bool("Color shading of interior of cells");
2390 if (prm.get("Color function") == "default")
2392 else if (prm.get("Color function") == "grey scale")
2394 else if (prm.get("Color function") == "reverse grey scale")
2396 else
2397 // we shouldn't get here, since the parameter object should already have
2398 // checked that the given value is valid
2399 Assert(false, ExcInternalError());
2400 }
2401
2402
2404 : compression_level(compression_level)
2405 {}
2406
2407
2408 TecplotFlags::TecplotFlags(const char *zone_name, const double solution_time)
2409 : zone_name(zone_name)
2410 , solution_time(solution_time)
2411 {}
2412
2413
2414
2415 std::size_t
2417 {
2418 return sizeof(*this) + MemoryConsumption::memory_consumption(zone_name);
2419 }
2420
2421
2422
2423 VtkFlags::VtkFlags(const double time,
2424 const unsigned int cycle,
2425 const bool print_date_and_time,
2426 const CompressionLevel compression_level,
2427 const bool write_higher_order_cells,
2428 const std::map<std::string, std::string> &physical_units)
2429 : time(time)
2430 , cycle(cycle)
2431 , print_date_and_time(print_date_and_time)
2432 , compression_level(compression_level)
2433 , write_higher_order_cells(write_higher_order_cells)
2434 , physical_units(physical_units)
2435 {}
2436
2437
2438
2440 parse_output_format(const std::string &format_name)
2441 {
2442 if (format_name == "none")
2443 return none;
2444
2445 if (format_name == "dx")
2446 return dx;
2447
2448 if (format_name == "ucd")
2449 return ucd;
2450
2451 if (format_name == "gnuplot")
2452 return gnuplot;
2453
2454 if (format_name == "povray")
2455 return povray;
2456
2457 if (format_name == "eps")
2458 return eps;
2459
2460 if (format_name == "gmv")
2461 return gmv;
2462
2463 if (format_name == "tecplot")
2464 return tecplot;
2465
2466 if (format_name == "vtk")
2467 return vtk;
2468
2469 if (format_name == "vtu")
2470 return vtu;
2471
2472 if (format_name == "deal.II intermediate")
2473 return deal_II_intermediate;
2474
2475 if (format_name == "hdf5")
2476 return hdf5;
2477
2478 AssertThrow(false,
2479 ExcMessage("The given file format name is not recognized: <" +
2480 format_name + ">"));
2481
2482 // return something invalid
2483 return OutputFormat(-1);
2484 }
2485
2486
2487
2488 std::string
2490 {
2491 return "none|dx|ucd|gnuplot|povray|eps|gmv|tecplot|vtk|vtu|hdf5|svg|deal.II intermediate";
2492 }
2493
2494
2495
2496 std::string
2497 default_suffix(const OutputFormat output_format)
2498 {
2499 switch (output_format)
2500 {
2501 case none:
2502 return "";
2503 case dx:
2504 return ".dx";
2505 case ucd:
2506 return ".inp";
2507 case gnuplot:
2508 return ".gnuplot";
2509 case povray:
2510 return ".pov";
2511 case eps:
2512 return ".eps";
2513 case gmv:
2514 return ".gmv";
2515 case tecplot:
2516 return ".dat";
2517 case vtk:
2518 return ".vtk";
2519 case vtu:
2520 return ".vtu";
2522 return ".d2";
2523 case hdf5:
2524 return ".h5";
2525 case svg:
2526 return ".svg";
2527 default:
2528 Assert(false, ExcNotImplemented());
2529 return "";
2530 }
2531 }
2532
2533
2534 //----------------------------------------------------------------------//
2535
2536
2541 template <int dim, int spacedim>
2542 std::vector<Point<spacedim>>
2543 get_node_positions(const std::vector<Patch<dim, spacedim>> &patches)
2544 {
2545 Assert(dim <= 3, ExcNotImplemented());
2546 static const std::array<unsigned int, 5> table = {{0, 1, 3, 2, 4}};
2547
2548 std::vector<Point<spacedim>> node_positions;
2549 for (const auto &patch : patches)
2550 {
2551 // special treatment of non-hypercube cells
2552 if (patch.reference_cell != ReferenceCells::get_hypercube<dim>())
2553 {
2554 for (unsigned int point_no = 0; point_no < patch.data.n_cols();
2555 ++point_no)
2556 node_positions.emplace_back(get_node_location(
2557 patch,
2559 table[point_no] :
2560 point_no)));
2561 }
2562 else
2563 {
2564 const unsigned int n_subdivisions = patch.n_subdivisions;
2565 const unsigned int n = n_subdivisions + 1;
2566
2567 switch (dim)
2568 {
2569 case 0:
2570 node_positions.emplace_back(
2571 get_equispaced_location(patch, {}, n_subdivisions));
2572 break;
2573 case 1:
2574 for (unsigned int i1 = 0; i1 < n; ++i1)
2575 node_positions.emplace_back(
2576 get_equispaced_location(patch, {i1}, n_subdivisions));
2577 break;
2578 case 2:
2579 for (unsigned int i2 = 0; i2 < n; ++i2)
2580 for (unsigned int i1 = 0; i1 < n; ++i1)
2581 node_positions.emplace_back(get_equispaced_location(
2582 patch, {i1, i2}, n_subdivisions));
2583 break;
2584 case 3:
2585 for (unsigned int i3 = 0; i3 < n; ++i3)
2586 for (unsigned int i2 = 0; i2 < n; ++i2)
2587 for (unsigned int i1 = 0; i1 < n; ++i1)
2588 node_positions.emplace_back(get_equispaced_location(
2589 patch, {i1, i2, i3}, n_subdivisions));
2590 break;
2591
2592 default:
2593 Assert(false, ExcInternalError());
2594 }
2595 }
2596 }
2597
2598 return node_positions;
2599 }
2600
2601
2602 template <int dim, int spacedim, typename StreamType>
2603 void
2604 write_nodes(const std::vector<Patch<dim, spacedim>> &patches, StreamType &out)
2605 {
2606 // Obtain the node locations, and then output them via the given stream
2607 // object
2608 const std::vector<Point<spacedim>> node_positions =
2609 get_node_positions(patches);
2610
2611 int count = 0;
2612 for (const auto &node : node_positions)
2613 out.write_point(count++, node);
2614 out.flush_points();
2615 }
2616
2617
2618
2619 template <int dim, int spacedim, typename StreamType>
2620 void
2621 write_cells(const std::vector<Patch<dim, spacedim>> &patches, StreamType &out)
2622 {
2623 Assert(dim <= 3, ExcNotImplemented());
2624 unsigned int count = 0;
2625 unsigned int first_vertex_of_patch = 0;
2626 for (const auto &patch : patches)
2627 {
2628 // special treatment of simplices since they are not subdivided
2629 if (patch.reference_cell != ReferenceCells::get_hypercube<dim>())
2630 {
2631 out.write_cell_single(count++,
2632 first_vertex_of_patch,
2633 patch.data.n_cols(),
2634 patch.reference_cell);
2635 first_vertex_of_patch += patch.data.n_cols();
2636 }
2637 else // hypercube cell
2638 {
2639 const unsigned int n_subdivisions = patch.n_subdivisions;
2640 const unsigned int n = n_subdivisions + 1;
2641
2642 switch (dim)
2643 {
2644 case 0:
2645 {
2646 const unsigned int offset = first_vertex_of_patch;
2647 out.template write_cell<0>(count++, offset, {});
2648 break;
2650
2651 case 1:
2652 {
2653 constexpr unsigned int d1 = 1;
2654
2655 for (unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
2656 {
2657 const unsigned int offset =
2658 first_vertex_of_patch + i1 * d1;
2659 out.template write_cell<1>(count++, offset, {{d1}});
2660 }
2661
2662 break;
2663 }
2664
2665 case 2:
2666 {
2667 constexpr unsigned int d1 = 1;
2668 const unsigned int d2 = n;
2670 for (unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
2671 for (unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
2672 {
2673 const unsigned int offset =
2674 first_vertex_of_patch + i2 * d2 + i1 * d1;
2675 out.template write_cell<2>(count++,
2676 offset,
2677 {{d1, d2}});
2678 }
2679
2680 break;
2681 }
2682
2683 case 3:
2684 {
2685 constexpr unsigned int d1 = 1;
2686 const unsigned int d2 = n;
2687 const unsigned int d3 = n * n;
2688
2689 for (unsigned int i3 = 0; i3 < n_subdivisions; ++i3)
2690 for (unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
2691 for (unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
2692 {
2693 const unsigned int offset = first_vertex_of_patch +
2694 i3 * d3 + i2 * d2 +
2695 i1 * d1;
2696 out.template write_cell<3>(count++,
2697 offset,
2698 {{d1, d2, d3}});
2699 }
2700
2701 break;
2702 }
2703 default:
2705 }
2706
2707 // Update the number of the first vertex of this patch
2708 first_vertex_of_patch +=
2709 Utilities::fixed_power<dim>(n_subdivisions + 1);
2710 }
2711 }
2712
2713 out.flush_cells();
2714 }
2715
2716
2717
2718 template <int dim, int spacedim, typename StreamType>
2719 void
2721 StreamType & out,
2722 const bool legacy_format)
2723 {
2724 Assert(dim <= 3 && dim > 1, ExcNotImplemented());
2725 unsigned int first_vertex_of_patch = 0;
2726 // Array to hold all the node numbers of a cell
2727 std::vector<unsigned> connectivity;
2728
2729 for (const auto &patch : patches)
2730 {
2731 if (patch.reference_cell != ReferenceCells::get_hypercube<dim>())
2732 {
2733 connectivity.resize(patch.data.n_cols());
2734
2735 for (unsigned int i = 0; i < patch.data.n_cols(); ++i)
2736 connectivity[i] = i;
2737
2738 out.template write_high_order_cell<dim>(first_vertex_of_patch,
2739 connectivity);
2740
2741 first_vertex_of_patch += patch.data.n_cols();
2742 }
2743 else
2744 {
2745 const unsigned int n_subdivisions = patch.n_subdivisions;
2746 const unsigned int n = n_subdivisions + 1;
2747
2748 connectivity.resize(Utilities::fixed_power<dim>(n));
2749
2750 switch (dim)
2751 {
2752 case 0:
2753 {
2754 Assert(false,
2755 ExcMessage("Point-like cells should not be possible "
2756 "when writing higher-order cells."));
2757 break;
2758 }
2759 case 1:
2760 {
2761 for (unsigned int i1 = 0; i1 < n_subdivisions + 1; ++i1)
2762 {
2763 const unsigned int local_index = i1;
2764 const unsigned int connectivity_index =
2765 patch.reference_cell
2766 .template vtk_lexicographic_to_node_index<1>(
2767 {{i1}}, {{n_subdivisions}}, legacy_format);
2768 connectivity[connectivity_index] = local_index;
2769 }
2770
2771 break;
2772 }
2773 case 2:
2774 {
2775 for (unsigned int i2 = 0; i2 < n_subdivisions + 1; ++i2)
2776 for (unsigned int i1 = 0; i1 < n_subdivisions + 1; ++i1)
2777 {
2778 const unsigned int local_index = i2 * n + i1;
2779 const unsigned int connectivity_index =
2780 patch.reference_cell
2781 .template vtk_lexicographic_to_node_index<2>(
2782 {{i1, i2}},
2783 {{n_subdivisions, n_subdivisions}},
2784 legacy_format);
2785 connectivity[connectivity_index] = local_index;
2786 }
2787
2788 break;
2789 }
2790 case 3:
2791 {
2792 for (unsigned int i3 = 0; i3 < n_subdivisions + 1; ++i3)
2793 for (unsigned int i2 = 0; i2 < n_subdivisions + 1; ++i2)
2794 for (unsigned int i1 = 0; i1 < n_subdivisions + 1; ++i1)
2795 {
2796 const unsigned int local_index =
2797 i3 * n * n + i2 * n + i1;
2798 const unsigned int connectivity_index =
2799 patch.reference_cell
2800 .template vtk_lexicographic_to_node_index<3>(
2801 {{i1, i2, i3}},
2802 {{n_subdivisions,
2803 n_subdivisions,
2804 n_subdivisions}},
2805 legacy_format);
2806 connectivity[connectivity_index] = local_index;
2807 }
2808
2809 break;
2810 }
2811 default:
2812 Assert(false, ExcNotImplemented());
2813 }
2814
2815 // Having so set up the 'connectivity' data structure,
2816 // output it:
2817 out.template write_high_order_cell<dim>(first_vertex_of_patch,
2818 connectivity);
2819
2820 // Finally update the number of the first vertex of this patch
2821 first_vertex_of_patch += Utilities::fixed_power<dim>(n);
2822 }
2823 }
2824
2825 out.flush_cells();
2826 }
2827
2828
2829 template <int dim, int spacedim, class StreamType>
2830 void
2831 write_data(const std::vector<Patch<dim, spacedim>> &patches,
2832 unsigned int n_data_sets,
2833 const bool double_precision,
2834 StreamType & out)
2835 {
2836 Assert(dim <= 3, ExcNotImplemented());
2837 unsigned int count = 0;
2838
2839 for (const auto &patch : patches)
2840 {
2841 const unsigned int n_subdivisions = patch.n_subdivisions;
2842 const unsigned int n = n_subdivisions + 1;
2843 // Length of loops in all dimensions
2844 Assert((patch.data.n_rows() == n_data_sets &&
2845 !patch.points_are_available) ||
2846 (patch.data.n_rows() == n_data_sets + spacedim &&
2847 patch.points_are_available),
2849 (n_data_sets + spacedim) :
2850 n_data_sets,
2851 patch.data.n_rows()));
2852 Assert(patch.data.n_cols() == Utilities::fixed_power<dim>(n),
2853 ExcInvalidDatasetSize(patch.data.n_cols(), n));
2854
2855 std::vector<float> floats(n_data_sets);
2856 std::vector<double> doubles(n_data_sets);
2858 // Data is already in lexicographic ordering
2859 for (unsigned int i = 0; i < Utilities::fixed_power<dim>(n);
2860 ++i, ++count)
2861 if (double_precision)
2862 {
2863 for (unsigned int data_set = 0; data_set < n_data_sets;
2864 ++data_set)
2865 doubles[data_set] = patch.data(data_set, i);
2866 out.write_dataset(count, doubles);
2867 }
2868 else
2869 {
2870 for (unsigned int data_set = 0; data_set < n_data_sets;
2871 ++data_set)
2872 floats[data_set] = patch.data(data_set, i);
2873 out.write_dataset(count, floats);
2874 }
2875 }
2876 }
2877
2878
2880 namespace
2881 {
2891 svg_project_point(Point<3> point,
2892 Point<3> camera_position,
2893 Point<3> camera_direction,
2894 Point<3> camera_horizontal,
2895 float camera_focus)
2896 {
2897 Point<3> camera_vertical;
2898 camera_vertical[0] = camera_horizontal[1] * camera_direction[2] -
2899 camera_horizontal[2] * camera_direction[1];
2900 camera_vertical[1] = camera_horizontal[2] * camera_direction[0] -
2901 camera_horizontal[0] * camera_direction[2];
2902 camera_vertical[2] = camera_horizontal[0] * camera_direction[1] -
2903 camera_horizontal[1] * camera_direction[0];
2904
2905 float phi;
2906 phi = camera_focus;
2907 phi /= (point[0] - camera_position[0]) * camera_direction[0] +
2908 (point[1] - camera_position[1]) * camera_direction[1] +
2909 (point[2] - camera_position[2]) * camera_direction[2];
2910
2911 Point<3> projection;
2912 projection[0] =
2913 camera_position[0] + phi * (point[0] - camera_position[0]);
2914 projection[1] =
2915 camera_position[1] + phi * (point[1] - camera_position[1]);
2916 projection[2] =
2917 camera_position[2] + phi * (point[2] - camera_position[2]);
2918
2919 Point<2> projection_decomposition;
2920 projection_decomposition[0] = (projection[0] - camera_position[0] -
2921 camera_focus * camera_direction[0]) *
2922 camera_horizontal[0];
2923 projection_decomposition[0] += (projection[1] - camera_position[1] -
2924 camera_focus * camera_direction[1]) *
2925 camera_horizontal[1];
2926 projection_decomposition[0] += (projection[2] - camera_position[2] -
2927 camera_focus * camera_direction[2]) *
2928 camera_horizontal[2];
2929
2930 projection_decomposition[1] = (projection[0] - camera_position[0] -
2931 camera_focus * camera_direction[0]) *
2932 camera_vertical[0];
2933 projection_decomposition[1] += (projection[1] - camera_position[1] -
2934 camera_focus * camera_direction[1]) *
2935 camera_vertical[1];
2936 projection_decomposition[1] += (projection[2] - camera_position[2] -
2937 camera_focus * camera_direction[2]) *
2938 camera_vertical[2];
2939
2940 return projection_decomposition;
2941 }
2942
2943
2948 Point<6>
2949 svg_get_gradient_parameters(Point<3> points[])
2950 {
2951 Point<3> v_min, v_max, v_inter;
2952
2953 // Use the Bubblesort algorithm to sort the points with respect to the
2954 // third coordinate
2955 for (int i = 0; i < 2; ++i)
2956 {
2957 for (int j = 0; j < 2 - i; ++j)
2958 {
2959 if (points[j][2] > points[j + 1][2])
2960 {
2961 Point<3> temp = points[j];
2962 points[j] = points[j + 1];
2963 points[j + 1] = temp;
2964 }
2965 }
2966 }
2967
2968 // save the related three-dimensional vectors v_min, v_inter, and v_max
2969 v_min = points[0];
2970 v_inter = points[1];
2971 v_max = points[2];
2972
2973 Point<2> A[2];
2975
2976 // determine the plane offset c
2977 A[0][0] = v_max[0] - v_min[0];
2978 A[0][1] = v_inter[0] - v_min[0];
2979 A[1][0] = v_max[1] - v_min[1];
2980 A[1][1] = v_inter[1] - v_min[1];
2981
2982 b[0] = -v_min[0];
2983 b[1] = -v_min[1];
2984
2985 double x, sum;
2986 bool col_change = false;
2987
2988 if (A[0][0] == 0)
2989 {
2990 col_change = true;
2991
2992 A[0][0] = A[0][1];
2993 A[0][1] = 0;
2994
2995 double temp = A[1][0];
2996 A[1][0] = A[1][1];
2997 A[1][1] = temp;
2998 }
2999
3000 for (unsigned int k = 0; k < 1; ++k)
3001 {
3002 for (unsigned int i = k + 1; i < 2; ++i)
3003 {
3004 x = A[i][k] / A[k][k];
3005
3006 for (unsigned int j = k + 1; j < 2; ++j)
3007 A[i][j] = A[i][j] - A[k][j] * x;
3008
3009 b[i] = b[i] - b[k] * x;
3010 }
3011 }
3012
3013 b[1] = b[1] / A[1][1];
3014
3015 for (int i = 0; i >= 0; i--)
3016 {
3017 sum = b[i];
3018
3019 for (unsigned int j = i + 1; j < 2; ++j)
3020 sum = sum - A[i][j] * b[j];
3021
3022 b[i] = sum / A[i][i];
3023 }
3024
3025 if (col_change)
3026 {
3027 double temp = b[0];
3028 b[0] = b[1];
3029 b[1] = temp;
3030 }
3031
3032 double c = b[0] * (v_max[2] - v_min[2]) + b[1] * (v_inter[2] - v_min[2]) +
3033 v_min[2];
3034
3035 // Determine the first entry of the gradient (phi, cf. documentation)
3036 A[0][0] = v_max[0] - v_min[0];
3037 A[0][1] = v_inter[0] - v_min[0];
3038 A[1][0] = v_max[1] - v_min[1];
3039 A[1][1] = v_inter[1] - v_min[1];
3040
3041 b[0] = 1.0 - v_min[0];
3042 b[1] = -v_min[1];
3043
3044 col_change = false;
3046 if (A[0][0] == 0)
3047 {
3048 col_change = true;
3049
3050 A[0][0] = A[0][1];
3051 A[0][1] = 0;
3052
3053 double temp = A[1][0];
3054 A[1][0] = A[1][1];
3055 A[1][1] = temp;
3056 }
3057
3058 for (unsigned int k = 0; k < 1; ++k)
3059 {
3060 for (unsigned int i = k + 1; i < 2; ++i)
3061 {
3062 x = A[i][k] / A[k][k];
3063
3064 for (unsigned int j = k + 1; j < 2; ++j)
3065 A[i][j] = A[i][j] - A[k][j] * x;
3066
3067 b[i] = b[i] - b[k] * x;
3068 }
3069 }
3070
3071 b[1] = b[1] / A[1][1];
3072
3073 for (int i = 0; i >= 0; i--)
3074 {
3075 sum = b[i];
3076
3077 for (unsigned int j = i + 1; j < 2; ++j)
3078 sum = sum - A[i][j] * b[j];
3079
3080 b[i] = sum / A[i][i];
3081 }
3082
3083 if (col_change)
3084 {
3085 double temp = b[0];
3086 b[0] = b[1];
3087 b[1] = temp;
3088 }
3089
3090 gradient[0] = b[0] * (v_max[2] - v_min[2]) +
3091 b[1] * (v_inter[2] - v_min[2]) - c + v_min[2];
3092
3093 // determine the second entry of the gradient
3094 A[0][0] = v_max[0] - v_min[0];
3095 A[0][1] = v_inter[0] - v_min[0];
3096 A[1][0] = v_max[1] - v_min[1];
3097 A[1][1] = v_inter[1] - v_min[1];
3098
3099 b[0] = -v_min[0];
3100 b[1] = 1.0 - v_min[1];
3101
3102 col_change = false;
3103
3104 if (A[0][0] == 0)
3106 col_change = true;
3107
3108 A[0][0] = A[0][1];
3109 A[0][1] = 0;
3110
3111 double temp = A[1][0];
3112 A[1][0] = A[1][1];
3113 A[1][1] = temp;
3114 }
3115
3116 for (unsigned int k = 0; k < 1; ++k)
3117 {
3118 for (unsigned int i = k + 1; i < 2; ++i)
3119 {
3120 x = A[i][k] / A[k][k];
3121
3122 for (unsigned int j = k + 1; j < 2; ++j)
3123 A[i][j] = A[i][j] - A[k][j] * x;
3124
3125 b[i] = b[i] - b[k] * x;
3126 }
3127 }
3128
3129 b[1] = b[1] / A[1][1];
3130
3131 for (int i = 0; i >= 0; i--)
3132 {
3133 sum = b[i];
3134
3135 for (unsigned int j = i + 1; j < 2; ++j)
3136 sum = sum - A[i][j] * b[j];
3137
3138 b[i] = sum / A[i][i];
3139 }
3140
3141 if (col_change)
3142 {
3143 double temp = b[0];
3144 b[0] = b[1];
3145 b[1] = temp;
3146 }
3147
3148 gradient[1] = b[0] * (v_max[2] - v_min[2]) +
3149 b[1] * (v_inter[2] - v_min[2]) - c + v_min[2];
3150
3151 // normalize the gradient
3152 gradient /= gradient.norm();
3153
3154 const double lambda = -gradient[0] * (v_min[0] - v_max[0]) -
3155 gradient[1] * (v_min[1] - v_max[1]);
3156
3157 Point<6> gradient_parameters;
3158
3159 gradient_parameters[0] = v_min[0];
3160 gradient_parameters[1] = v_min[1];
3161
3162 gradient_parameters[2] = v_min[0] + lambda * gradient[0];
3163 gradient_parameters[3] = v_min[1] + lambda * gradient[1];
3164
3165 gradient_parameters[4] = v_min[2];
3166 gradient_parameters[5] = v_max[2];
3167
3168 return gradient_parameters;
3169 }
3170 } // namespace
3171
3172
3173
3174 template <int dim, int spacedim>
3175 void
3177 const std::vector<Patch<dim, spacedim>> &patches,
3178 const std::vector<std::string> & data_names,
3179 const std::vector<
3180 std::tuple<unsigned int,
3181 unsigned int,
3182 std::string,
3184 const UcdFlags &flags,
3185 std::ostream & out)
3186 {
3187 // Note that while in theory dim==0 should be implemented, this is not
3188 // tested, therefore currently not allowed.
3189 AssertThrow(dim > 0, ExcNotImplemented());
3190
3191 AssertThrow(out.fail() == false, ExcIO());
3192
3193#ifndef DEAL_II_WITH_MPI
3194 // verify that there are indeed patches to be written out. most of the
3195 // times, people just forget to call build_patches when there are no
3196 // patches, so a warning is in order. that said, the assertion is disabled
3197 // if we support MPI since then it can happen that on the coarsest mesh, a
3198 // processor simply has no cells it actually owns, and in that case it is
3199 // legit if there are no patches
3200 Assert(patches.size() > 0, ExcNoPatches());
3201#else
3202 if (patches.size() == 0)
3203 return;
3204#endif
3205
3206 const unsigned int n_data_sets = data_names.size();
3207
3208 UcdStream ucd_out(out, flags);
3209
3210 // first count the number of cells and cells for later use
3211 unsigned int n_nodes;
3212 unsigned int n_cells;
3213 std::tie(n_nodes, n_cells) = count_nodes_and_cells(patches);
3214 //---------------------
3215 // preamble
3216 if (flags.write_preamble)
3217 {
3218 out
3219 << "# This file was generated by the deal.II library." << '\n'
3220 << "# Date = " << Utilities::System::get_date() << '\n'
3221 << "# Time = " << Utilities::System::get_time() << '\n'
3222 << "#" << '\n'
3223 << "# For a description of the UCD format see the AVS Developer's guide."
3224 << '\n'
3225 << "#" << '\n';
3226 }
3227
3228 // start with ucd data
3229 out << n_nodes << ' ' << n_cells << ' ' << n_data_sets << ' ' << 0
3230 << ' ' // no cell data at present
3231 << 0 // no model data
3232 << '\n';
3233
3234 write_nodes(patches, ucd_out);
3235 out << '\n';
3236
3237 write_cells(patches, ucd_out);
3238 out << '\n';
3239
3240 //---------------------------
3241 // now write data
3242 if (n_data_sets != 0)
3243 {
3244 out << n_data_sets << " "; // number of vectors
3245 for (unsigned int i = 0; i < n_data_sets; ++i)
3246 out << 1 << ' '; // number of components;
3247 // only 1 supported presently
3248 out << '\n';
3249
3250 for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
3251 out << data_names[data_set]
3252 << ",dimensionless" // no units supported at present
3253 << '\n';
3254
3255 write_data(patches, n_data_sets, true, ucd_out);
3256 }
3257 // make sure everything now gets to disk
3258 out.flush();
3259
3260 // assert the stream is still ok
3261 AssertThrow(out.fail() == false, ExcIO());
3262 }
3263
3264
3265 template <int dim, int spacedim>
3266 void
3268 const std::vector<Patch<dim, spacedim>> &patches,
3269 const std::vector<std::string> & data_names,
3270 const std::vector<
3271 std::tuple<unsigned int,
3272 unsigned int,
3273 std::string,
3275 const DXFlags &flags,
3276 std::ostream & out)
3277 {
3278 // Point output is currently not implemented.
3279 AssertThrow(dim > 0, ExcNotImplemented());
3280
3281 AssertThrow(out.fail() == false, ExcIO());
3282
3283#ifndef DEAL_II_WITH_MPI
3284 // verify that there are indeed patches to be written out. most of the
3285 // times, people just forget to call build_patches when there are no
3286 // patches, so a warning is in order. that said, the assertion is disabled
3287 // if we support MPI since then it can happen that on the coarsest mesh, a
3288 // processor simply has no cells it actually owns, and in that case it is
3289 // legit if there are no patches
3290 Assert(patches.size() > 0, ExcNoPatches());
3291#else
3292 if (patches.size() == 0)
3293 return;
3294#endif
3295 // Stream with special features for dx output
3296 DXStream dx_out(out, flags);
3297
3298 // Variable counting the offset of binary data.
3299 unsigned int offset = 0;
3300
3301 const unsigned int n_data_sets = data_names.size();
3302
3303 // first count the number of cells and cells for later use
3304 unsigned int n_nodes;
3305 unsigned int n_cells;
3306 std::tie(n_nodes, n_cells) = count_nodes_and_cells(patches);
3307
3308 // start with vertices order is lexicographical, x varying fastest
3309 out << "object \"vertices\" class array type float rank 1 shape "
3310 << spacedim << " items " << n_nodes;
3311
3312 if (flags.coordinates_binary)
3313 {
3314 out << " lsb ieee data 0" << '\n';
3315 offset += n_nodes * spacedim * sizeof(float);
3316 }
3317 else
3318 {
3319 out << " data follows" << '\n';
3320 write_nodes(patches, dx_out);
3321 }
3322
3323 //-----------------------------
3324 // first write the coordinates of all vertices
3325
3326 //---------------------------------------
3327 // write cells
3328 out << "object \"cells\" class array type int rank 1 shape "
3329 << GeometryInfo<dim>::vertices_per_cell << " items " << n_cells;
3330
3331 if (flags.int_binary)
3332 {
3333 out << " lsb binary data " << offset << '\n';
3334 offset += n_cells * sizeof(int);
3335 }
3336 else
3337 {
3338 out << " data follows" << '\n';
3339 write_cells(patches, dx_out);
3340 out << '\n';
3341 }
3342
3343
3344 out << "attribute \"element type\" string \"";
3345 if (dim == 1)
3346 out << "lines";
3347 if (dim == 2)
3348 out << "quads";
3349 if (dim == 3)
3350 out << "cubes";
3351 out << "\"" << '\n' << "attribute \"ref\" string \"positions\"" << '\n';
3352
3353 // TODO:[GK] Patches must be of same size!
3354 //---------------------------
3355 // write neighbor information
3356 if (flags.write_neighbors)
3357 {
3358 out << "object \"neighbors\" class array type int rank 1 shape "
3359 << GeometryInfo<dim>::faces_per_cell << " items " << n_cells
3360 << " data follows";
3361
3362 for (const auto &patch : patches)
3363 {
3364 const unsigned int n = patch.n_subdivisions;
3365 const unsigned int n1 = (dim > 0) ? n : 1;
3366 const unsigned int n2 = (dim > 1) ? n : 1;
3367 const unsigned int n3 = (dim > 2) ? n : 1;
3368 const unsigned int x_minus = (dim > 0) ? 0 : 0;
3369 const unsigned int x_plus = (dim > 0) ? 1 : 0;
3370 const unsigned int y_minus = (dim > 1) ? 2 : 0;
3371 const unsigned int y_plus = (dim > 1) ? 3 : 0;
3372 const unsigned int z_minus = (dim > 2) ? 4 : 0;
3373 const unsigned int z_plus = (dim > 2) ? 5 : 0;
3374 unsigned int cells_per_patch = Utilities::fixed_power<dim>(n);
3375 unsigned int dx = 1;
3376 unsigned int dy = n;
3377 unsigned int dz = n * n;
3378
3379 const unsigned int patch_start =
3380 patch.patch_index * cells_per_patch;
3381
3382 for (unsigned int i3 = 0; i3 < n3; ++i3)
3383 for (unsigned int i2 = 0; i2 < n2; ++i2)
3384 for (unsigned int i1 = 0; i1 < n1; ++i1)
3385 {
3386 const unsigned int nx = i1 * dx;
3387 const unsigned int ny = i2 * dy;
3388 const unsigned int nz = i3 * dz;
3389
3390 // There are no neighbors for dim==0. Note that this case is
3391 // caught by the AssertThrow at the beginning of this
3392 // function anyway. This condition avoids compiler warnings.
3393 if (dim < 1)
3394 continue;
3395
3396 out << '\n';
3397 // Direction -x Last cell in row of other patch
3398 if (i1 == 0)
3399 {
3400 const unsigned int nn = patch.neighbors[x_minus];
3401 out << '\t';
3402 if (nn != patch.no_neighbor)
3403 out
3404 << (nn * cells_per_patch + ny + nz + dx * (n - 1));
3405 else
3406 out << "-1";
3407 }
3408 else
3409 {
3410 out << '\t' << patch_start + nx - dx + ny + nz;
3411 }
3412 // Direction +x First cell in row of other patch
3413 if (i1 == n - 1)
3414 {
3415 const unsigned int nn = patch.neighbors[x_plus];
3416 out << '\t';
3417 if (nn != patch.no_neighbor)
3418 out << (nn * cells_per_patch + ny + nz);
3419 else
3420 out << "-1";
3421 }
3422 else
3423 {
3424 out << '\t' << patch_start + nx + dx + ny + nz;
3425 }
3426 if (dim < 2)
3427 continue;
3428 // Direction -y
3429 if (i2 == 0)
3430 {
3431 const unsigned int nn = patch.neighbors[y_minus];
3432 out << '\t';
3433 if (nn != patch.no_neighbor)
3434 out
3435 << (nn * cells_per_patch + nx + nz + dy * (n - 1));
3436 else
3437 out << "-1";
3438 }
3439 else
3440 {
3441 out << '\t' << patch_start + nx + ny - dy + nz;
3442 }
3443 // Direction +y
3444 if (i2 == n - 1)
3445 {
3446 const unsigned int nn = patch.neighbors[y_plus];
3447 out << '\t';
3448 if (nn != patch.no_neighbor)
3449 out << (nn * cells_per_patch + nx + nz);
3450 else
3451 out << "-1";
3452 }
3453 else
3454 {
3455 out << '\t' << patch_start + nx + ny + dy + nz;
3456 }
3457 if (dim < 3)
3458 continue;
3459
3460 // Direction -z
3461 if (i3 == 0)
3462 {
3463 const unsigned int nn = patch.neighbors[z_minus];
3464 out << '\t';
3465 if (nn != patch.no_neighbor)
3466 out
3467 << (nn * cells_per_patch + nx + ny + dz * (n - 1));
3468 else
3469 out << "-1";
3470 }
3471 else
3472 {
3473 out << '\t' << patch_start + nx + ny + nz - dz;
3474 }
3475 // Direction +z
3476 if (i3 == n - 1)
3477 {
3478 const unsigned int nn = patch.neighbors[z_plus];
3479 out << '\t';
3480 if (nn != patch.no_neighbor)
3481 out << (nn * cells_per_patch + nx + ny);
3482 else
3483 out << "-1";
3484 }
3485 else
3486 {
3487 out << '\t' << patch_start + nx + ny + nz + dz;
3488 }
3489 }
3490 out << '\n';
3491 }
3492 }
3493 //---------------------------
3494 // now write data
3495 if (n_data_sets != 0)
3496 {
3497 out << "object \"data\" class array type float rank 1 shape "
3498 << n_data_sets << " items " << n_nodes;
3499
3500 if (flags.data_binary)
3501 {
3502 out << " lsb ieee data " << offset << '\n';
3503 offset += n_data_sets * n_nodes *
3504 ((flags.data_double) ? sizeof(double) : sizeof(float));
3505 }
3506 else
3507 {
3508 out << " data follows" << '\n';
3509 write_data(patches, n_data_sets, flags.data_double, dx_out);
3510 }
3511
3512 // loop over all patches
3513 out << "attribute \"dep\" string \"positions\"" << '\n';
3514 }
3515 else
3516 {
3517 out << "object \"data\" class constantarray type float rank 0 items "
3518 << n_nodes << " data follows" << '\n'
3519 << '0' << '\n';
3520 }
3521
3522 // no model data
3523
3524 out << "object \"deal data\" class field" << '\n'
3525 << "component \"positions\" value \"vertices\"" << '\n'
3526 << "component \"connections\" value \"cells\"" << '\n'
3527 << "component \"data\" value \"data\"" << '\n';
3528
3529 if (flags.write_neighbors)
3530 out << "component \"neighbors\" value \"neighbors\"" << '\n';
3531
3532 {
3533 out << "attribute \"created\" string \"" << Utilities::System::get_date()
3534 << ' ' << Utilities::System::get_time() << '"' << '\n';
3535 }
3536
3537 out << "end" << '\n';
3538 // Write all binary data now
3539 if (flags.coordinates_binary)
3540 write_nodes(patches, dx_out);
3541 if (flags.int_binary)
3542 write_cells(patches, dx_out);
3543 if (flags.data_binary)
3544 write_data(patches, n_data_sets, flags.data_double, dx_out);
3545
3546 // make sure everything now gets to disk
3547 out.flush();
3548
3549 // assert the stream is still ok
3550 AssertThrow(out.fail() == false, ExcIO());
3551 }
3552
3553
3554
3555 template <int dim, int spacedim>
3556 void
3558 const std::vector<Patch<dim, spacedim>> &patches,
3559 const std::vector<std::string> & data_names,
3560 const std::vector<
3561 std::tuple<unsigned int,
3562 unsigned int,
3563 std::string,
3565 const GnuplotFlags &flags,
3566 std::ostream & out)
3567 {
3568 AssertThrow(out.fail() == false, ExcIO());
3569
3570#ifndef DEAL_II_WITH_MPI
3571 // verify that there are indeed patches to be written out. most
3572 // of the times, people just forget to call build_patches when there
3573 // are no patches, so a warning is in order. that said, the
3574 // assertion is disabled if we support MPI since then it can
3575 // happen that on the coarsest mesh, a processor simply has no
3576 // cells it actually owns, and in that case it is legit if there
3577 // are no patches
3578 Assert(patches.size() > 0, ExcNoPatches());
3579#else
3580 if (patches.size() == 0)
3581 return;
3582#endif
3583
3584 const unsigned int n_data_sets = data_names.size();
3585
3586 // write preamble
3587 {
3588 out << "# This file was generated by the deal.II library." << '\n'
3589 << "# Date = " << Utilities::System::get_date() << '\n'
3590 << "# Time = " << Utilities::System::get_time() << '\n'
3591 << "#" << '\n'
3592 << "# For a description of the GNUPLOT format see the GNUPLOT manual."
3593 << '\n'
3594 << "#" << '\n'
3595 << "# ";
3596
3597 AssertThrow(spacedim <= flags.space_dimension_labels.size(),
3599 for (unsigned int spacedim_n = 0; spacedim_n < spacedim; ++spacedim_n)
3600 {
3601 out << '<' << flags.space_dimension_labels.at(spacedim_n) << "> ";
3602 }
3603
3604 for (const auto &data_name : data_names)
3605 out << '<' << data_name << "> ";
3606 out << '\n';
3607 }
3608
3609
3610 // loop over all patches
3611 for (const auto &patch : patches)
3612 {
3613 const unsigned int n_subdivisions = patch.n_subdivisions;
3614 const unsigned int n_points_per_direction = n_subdivisions + 1;
3615
3616 Assert((patch.data.n_rows() == n_data_sets &&
3617 !patch.points_are_available) ||
3618 (patch.data.n_rows() == n_data_sets + spacedim &&
3619 patch.points_are_available),
3621 (n_data_sets + spacedim) :
3622 n_data_sets,
3623 patch.data.n_rows()));
3624
3625 auto output_point_data =
3626 [&out, &patch, n_data_sets](const unsigned int point_index) mutable {
3627 for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
3628 out << patch.data(data_set, point_index) << ' ';
3629 };
3630
3631 switch (dim)
3632 {
3633 case 0:
3634 {
3635 Assert(patch.reference_cell == ReferenceCells::Vertex,
3637 Assert(patch.data.n_cols() == 1,
3638 ExcInvalidDatasetSize(patch.data.n_cols(),
3639 n_subdivisions + 1));
3640
3641
3642 // compute coordinates for this patch point
3643 out << get_equispaced_location(patch, {}, n_subdivisions)
3644 << ' ';
3645 output_point_data(0);
3646 out << '\n';
3647 out << '\n';
3648 break;
3649 }
3650
3651 case 1:
3652 {
3653 Assert(patch.reference_cell == ReferenceCells::Line,
3655 Assert(patch.data.n_cols() ==
3656 Utilities::fixed_power<dim>(n_points_per_direction),
3657 ExcInvalidDatasetSize(patch.data.n_cols(),
3658 n_subdivisions + 1));
3659
3660 for (unsigned int i1 = 0; i1 < n_points_per_direction; ++i1)
3661 {
3662 // compute coordinates for this patch point
3663 out << get_equispaced_location(patch, {i1}, n_subdivisions)
3664 << ' ';
3665
3666 output_point_data(i1);
3667 out << '\n';
3668 }
3669 // end of patch
3670 out << '\n';
3671 out << '\n';
3672 break;
3673 }
3674
3675 case 2:
3676 {
3677 if (patch.reference_cell == ReferenceCells::Quadrilateral)
3678 {
3679 Assert(patch.data.n_cols() == Utilities::fixed_power<dim>(
3680 n_points_per_direction),
3681 ExcInvalidDatasetSize(patch.data.n_cols(),
3682 n_subdivisions + 1));
3683
3684 for (unsigned int i2 = 0; i2 < n_points_per_direction; ++i2)
3685 {
3686 for (unsigned int i1 = 0; i1 < n_points_per_direction;
3687 ++i1)
3688 {
3689 // compute coordinates for this patch point
3690 out << get_equispaced_location(patch,
3691 {i1, i2},
3692 n_subdivisions)
3693 << ' ';
3694
3695 output_point_data(i1 + i2 * n_points_per_direction);
3696 out << '\n';
3697 }
3698 // end of row in patch
3699 out << '\n';
3700 }
3701 }
3702 else if (patch.reference_cell == ReferenceCells::Triangle)
3703 {
3704 Assert(n_subdivisions == 1, ExcNotImplemented());
3705
3706 Assert(patch.data.n_cols() == 3, ExcInternalError());
3707
3708 // Gnuplot can only plot surfaces if each facet of the
3709 // surface is a bilinear patch, or a subdivided bilinear
3710 // patch with equally many points along each row of the
3711 // subdivision. This is what the code above for
3712 // quadrilaterals does. We emulate this by repeating the
3713 // third point of a triangle twice so that there are two
3714 // points for that row as well -- i.e., we write a 2x2
3715 // bilinear patch where two of the points are collapsed onto
3716 // one vertex.
3717 //
3718 // This also matches the example here:
3719 // https://stackoverflow.com/questions/42784369/drawing-triangular-mesh-using-gnuplot
3720 out << get_node_location(patch, 0) << ' ';
3721 output_point_data(0);
3722 out << '\n';
3723
3724 out << get_node_location(patch, 1) << ' ';
3725 output_point_data(1);
3726 out << '\n';
3727 out << '\n'; // end of one row of points
3728
3729 out << get_node_location(patch, 2) << ' ';
3730 output_point_data(2);
3731 out << '\n';
3732
3733 out << get_node_location(patch, 2) << ' ';
3734 output_point_data(2);
3735 out << '\n';
3736 out << '\n'; // end of the second row of points
3737 out << '\n'; // end of the entire patch
3738 }
3739 else
3740 // There aren't any other reference cells in 2d than the
3741 // quadrilateral and the triangle. So whatever we got here
3742 // can't be any good
3743 Assert(false, ExcInternalError());
3744 // end of patch
3745 out << '\n';
3746
3747 break;
3748 }
3749
3750 case 3:
3751 {
3752 if (patch.reference_cell == ReferenceCells::Hexahedron)
3753 {
3754 Assert(patch.data.n_cols() == Utilities::fixed_power<dim>(
3755 n_points_per_direction),
3756 ExcInvalidDatasetSize(patch.data.n_cols(),
3757 n_subdivisions + 1));
3758
3759 // for all grid points: draw lines into all positive
3760 // coordinate directions if there is another grid point
3761 // there
3762 for (unsigned int i3 = 0; i3 < n_points_per_direction; ++i3)
3763 for (unsigned int i2 = 0; i2 < n_points_per_direction;
3764 ++i2)
3765 for (unsigned int i1 = 0; i1 < n_points_per_direction;
3766 ++i1)
3767 {
3768 // compute coordinates for this patch point
3769 const Point<spacedim> this_point =
3770 get_equispaced_location(patch,
3771 {i1, i2, i3},
3772 n_subdivisions);
3773 // line into positive x-direction if possible
3774 if (i1 < n_subdivisions)
3775 {
3776 // write point here and its data
3777 out << this_point << ' ';
3778 output_point_data(i1 +
3779 i2 * n_points_per_direction +
3780 i3 * n_points_per_direction *
3781 n_points_per_direction);
3782 out << '\n';
3783
3784 // write point there and its data
3785 out << get_equispaced_location(patch,
3786 {i1 + 1, i2, i3},
3787 n_subdivisions)
3788 << ' ';
3789
3790 output_point_data((i1 + 1) +
3791 i2 * n_points_per_direction +
3792 i3 * n_points_per_direction *
3793 n_points_per_direction);
3794 out << '\n';
3795
3796 // end of line
3797 out << '\n' << '\n';
3798 }
3799
3800 // line into positive y-direction if possible
3801 if (i2 < n_subdivisions)
3802 {
3803 // write point here and its data
3804 out << this_point << ' ';
3805 output_point_data(i1 +
3806 i2 * n_points_per_direction +
3807 i3 * n_points_per_direction *
3808 n_points_per_direction);
3809 out << '\n';
3810
3811 // write point there and its data
3812 out << get_equispaced_location(patch,
3813 {i1, i2 + 1, i3},
3814 n_subdivisions)
3815 << ' ';
3816
3817 output_point_data(
3818 i1 + (i2 + 1) * n_points_per_direction +
3819 i3 * n_points_per_direction *
3820 n_points_per_direction);
3821 out << '\n';
3822
3823 // end of line
3824 out << '\n' << '\n';
3825 }
3826
3827 // line into positive z-direction if possible
3828 if (i3 < n_subdivisions)
3829 {
3830 // write point here and its data
3831 out << this_point << ' ';
3832 output_point_data(i1 +
3833 i2 * n_points_per_direction +
3834 i3 * n_points_per_direction *
3835 n_points_per_direction);
3836 out << '\n';
3837
3838 // write point there and its data
3839 out << get_equispaced_location(patch,
3840 {i1, i2, i3 + 1},
3841 n_subdivisions)
3842 << ' ';
3843
3844 output_point_data(
3845 i1 + i2 * n_points_per_direction +
3846 (i3 + 1) * n_points_per_direction *
3847 n_points_per_direction);
3848 out << '\n';
3849 // end of line
3850 out << '\n' << '\n';
3851 }
3852 }
3853 }
3854 else if (patch.reference_cell == ReferenceCells::Tetrahedron)
3855 {
3856 Assert(n_subdivisions == 1, ExcNotImplemented());
3857
3858 // Draw the tetrahedron as a collection of two lines.
3859 for (const unsigned int v : {0, 1, 2, 0, 3, 2})
3860 {
3861 out << get_node_location(patch, v) << ' ';
3862 output_point_data(v);
3863 out << '\n';
3864 }
3865 out << '\n'; // end of first line
3866
3867 for (const unsigned int v : {3, 1})
3868 {
3869 out << get_node_location(patch, v) << ' ';
3870 output_point_data(v);
3871 out << '\n';
3872 }
3873 out << '\n'; // end of second line
3874 }
3875 else if (patch.reference_cell == ReferenceCells::Pyramid)
3876 {
3877 Assert(n_subdivisions == 1, ExcNotImplemented());
3878
3879 // Draw the pyramid as a collection of two lines.
3880 for (const unsigned int v : {0, 1, 3, 2, 0, 4, 1})
3881 {
3882 out << get_node_location(patch, v) << ' ';
3883 output_point_data(v);
3884 out << '\n';
3885 }
3886 out << '\n'; // end of first line
3887
3888 for (const unsigned int v : {2, 4, 3})
3889 {
3890 out << get_node_location(patch, v) << ' ';
3891 output_point_data(v);
3892 out << '\n';
3893 }
3894 out << '\n'; // end of second line
3895 }
3896 else if (patch.reference_cell == ReferenceCells::Wedge)
3897 {
3898 Assert(n_subdivisions == 1, ExcNotImplemented());
3899
3900 // Draw the wedge as a collection of three
3901 // lines. The first one wraps around the base,
3902 // goes up to the top, and wraps around that. The
3903 // second and third are just individual lines
3904 // going from base to top.
3905 for (const unsigned int v : {0, 1, 2, 0, 3, 4, 5, 3})
3906 {
3907 out << get_node_location(patch, v) << ' ';
3908 output_point_data(v);
3909 out << '\n';
3910 }
3911 out << '\n'; // end of first line
3912
3913 for (const unsigned int v : {1, 4})
3914 {
3915 out << get_node_location(patch, v) << ' ';
3916 output_point_data(v);
3917 out << '\n';
3918 }
3919 out << '\n'; // end of second line
3920
3921 for (const unsigned int v : {2, 5})
3922 {
3923 out << get_node_location(patch, v) << ' ';
3924 output_point_data(v);
3925 out << '\n';
3926 }
3927 out << '\n'; // end of second line
3928 }
3929 else
3930 // No other reference cells are currently implemented
3931 Assert(false, ExcNotImplemented());
3932
3933 break;
3934 }
3935
3936 default:
3937 Assert(false, ExcNotImplemented());
3938 }
3939 }
3940 // make sure everything now gets to disk
3941 out.flush();
3942
3943 AssertThrow(out.fail() == false, ExcIO());
3944 }
3945
3946
3947 namespace
3948 {
3949 template <int dim, int spacedim>
3950 void
3951 do_write_povray(const std::vector<Patch<dim, spacedim>> &,
3952 const std::vector<std::string> &,
3953 const PovrayFlags &,
3954 std::ostream &)
3955 {
3956 Assert(false,
3957 ExcMessage("Writing files in POVRAY format is only supported "
3958 "for two-dimensional meshes."));
3959 }
3960
3961
3962
3963 void
3964 do_write_povray(const std::vector<Patch<2, 2>> &patches,
3965 const std::vector<std::string> &data_names,
3966 const PovrayFlags & flags,
3967 std::ostream & out)
3968 {
3969 AssertThrow(out.fail() == false, ExcIO());
3970
3971#ifndef DEAL_II_WITH_MPI
3972 // verify that there are indeed patches to be written out. most
3973 // of the times, people just forget to call build_patches when there
3974 // are no patches, so a warning is in order. that said, the
3975 // assertion is disabled if we support MPI since then it can
3976 // happen that on the coarsest mesh, a processor simply has no cells it
3977 // actually owns, and in that case it is legit if there are no patches
3978 Assert(patches.size() > 0, ExcNoPatches());
3979#else
3980 if (patches.size() == 0)
3981 return;
3982#endif
3983 constexpr int dim = 2;
3984 (void)dim;
3985 constexpr int spacedim = 2;
3986
3987 const unsigned int n_data_sets = data_names.size();
3988 (void)n_data_sets;
3989
3990 // write preamble
3991 {
3992 out
3993 << "/* This file was generated by the deal.II library." << '\n'
3994 << " Date = " << Utilities::System::get_date() << '\n'
3995 << " Time = " << Utilities::System::get_time() << '\n'
3996 << '\n'
3997 << " For a description of the POVRAY format see the POVRAY manual."
3998 << '\n'
3999 << "*/ " << '\n';
4000
4001 // include files
4002 out << "#include \"colors.inc\" " << '\n'
4003 << "#include \"textures.inc\" " << '\n';
4004
4005
4006 // use external include file for textures, camera and light
4007 if (flags.external_data)
4008 out << "#include \"data.inc\" " << '\n';
4009 else // all definitions in data file
4010 {
4011 // camera
4012 out << '\n'
4013 << '\n'
4014 << "camera {" << '\n'
4015 << " location <1,4,-7>" << '\n'
4016 << " look_at <0,0,0>" << '\n'
4017 << " angle 30" << '\n'
4018 << "}" << '\n';
4019
4020 // light
4021 out << '\n'
4022 << "light_source {" << '\n'
4023 << " <1,4,-7>" << '\n'
4024 << " color Grey" << '\n'
4025 << "}" << '\n';
4026 out << '\n'
4027 << "light_source {" << '\n'
4028 << " <0,20,0>" << '\n'
4029 << " color White" << '\n'
4030 << "}" << '\n';
4031 }
4032 }
4033
4034 // max. and min. height of solution
4035 Assert(patches.size() > 0, ExcNoPatches());
4036 double hmin = patches[0].data(0, 0);
4037 double hmax = patches[0].data(0, 0);
4038
4039 for (const auto &patch : patches)
4040 {
4041 const unsigned int n_subdivisions = patch.n_subdivisions;
4042
4043 Assert((patch.data.n_rows() == n_data_sets &&
4044 !patch.points_are_available) ||
4045 (patch.data.n_rows() == n_data_sets + spacedim &&
4046 patch.points_are_available),
4048 (n_data_sets + spacedim) :
4049 n_data_sets,
4050 patch.data.n_rows()));
4051 Assert(patch.data.n_cols() ==
4052 Utilities::fixed_power<dim>(n_subdivisions + 1),
4053 ExcInvalidDatasetSize(patch.data.n_cols(),
4054 n_subdivisions + 1));
4055
4056 for (unsigned int i = 0; i < n_subdivisions + 1; ++i)
4057 for (unsigned int j = 0; j < n_subdivisions + 1; ++j)
4058 {
4059 const int dl = i * (n_subdivisions + 1) + j;
4060 if (patch.data(0, dl) < hmin)
4061 hmin = patch.data(0, dl);
4062 if (patch.data(0, dl) > hmax)
4063 hmax = patch.data(0, dl);
4064 }
4065 }
4066
4067 out << "#declare HMIN=" << hmin << ";" << '\n'
4068 << "#declare HMAX=" << hmax << ";" << '\n'
4069 << '\n';
4070
4071 if (!flags.external_data)
4072 {
4073 // texture with scaled niveau lines 10 lines in the surface
4074 out << "#declare Tex=texture{" << '\n'
4075 << " pigment {" << '\n'
4076 << " gradient y" << '\n'
4077 << " scale y*(HMAX-HMIN)*" << 0.1 << '\n'
4078 << " color_map {" << '\n'
4079 << " [0.00 color Light_Purple] " << '\n'
4080 << " [0.95 color Light_Purple] " << '\n'
4081 << " [1.00 color White] " << '\n'
4082 << "} } }" << '\n'
4083 << '\n';
4084 }
4085
4086 if (!flags.bicubic_patch)
4087 {
4088 // start of mesh header
4089 out << '\n' << "mesh {" << '\n';
4090 }
4091
4092 // loop over all patches
4093 for (const auto &patch : patches)
4094 {
4095 const unsigned int n_subdivisions = patch.n_subdivisions;
4096 const unsigned int n = n_subdivisions + 1;
4097 const unsigned int d1 = 1;
4098 const unsigned int d2 = n;
4099
4100 Assert((patch.data.n_rows() == n_data_sets &&
4101 !patch.points_are_available) ||
4102 (patch.data.n_rows() == n_data_sets + spacedim &&
4103 patch.points_are_available),
4105 (n_data_sets + spacedim) :
4106 n_data_sets,
4107 patch.data.n_rows()));
4108 Assert(patch.data.n_cols() == Utilities::fixed_power<dim>(n),
4109 ExcInvalidDatasetSize(patch.data.n_cols(),
4110 n_subdivisions + 1));
4111
4112
4113 std::vector<Point<spacedim>> ver(n * n);
4114
4115 for (unsigned int i2 = 0; i2 < n; ++i2)
4116 for (unsigned int i1 = 0; i1 < n; ++i1)
4117 {
4118 // compute coordinates for this patch point, storing in ver
4119 ver[i1 * d1 + i2 * d2] =
4120 get_equispaced_location(patch, {i1, i2}, n_subdivisions);
4121 }
4122
4123
4124 if (!flags.bicubic_patch)
4125 {
4126 // approximate normal vectors in patch
4127 std::vector<Point<3>> nrml;
4128 // only if smooth triangles are used
4129 if (flags.smooth)
4130 {
4131 nrml.resize(n * n);
4132 // These are difference quotients of the surface
4133 // mapping. We take them symmetric inside the
4134 // patch and one-sided at the edges
4135 Point<3> h1, h2;
4136 // Now compute normals in every point
4137 for (unsigned int i = 0; i < n; ++i)
4138 for (unsigned int j = 0; j < n; ++j)
4139 {
4140 const unsigned int il = (i == 0) ? i : (i - 1);
4141 const unsigned int ir =
4142 (i == n_subdivisions) ? i : (i + 1);
4143 const unsigned int jl = (j == 0) ? j : (j - 1);
4144 const unsigned int jr =
4145 (j == n_subdivisions) ? j : (j + 1);
4146
4147 h1(0) =
4148 ver[ir * d1 + j * d2](0) - ver[il * d1 + j * d2](0);
4149 h1(1) = patch.data(0, ir * d1 + j * d2) -
4150 patch.data(0, il * d1 + j * d2);
4151 h1(2) =
4152 ver[ir * d1 + j * d2](1) - ver[il * d1 + j * d2](1);
4153
4154 h2(0) =
4155 ver[i * d1 + jr * d2](0) - ver[i * d1 + jl * d2](0);
4156 h2(1) = patch.data(0, i * d1 + jr * d2) -
4157 patch.data(0, i * d1 + jl * d2);
4158 h2(2) =
4159 ver[i * d1 + jr * d2](1) - ver[i * d1 + jl * d2](1);
4160
4161 nrml[i * d1 + j * d2](0) =
4162 h1(1) * h2(2) - h1(2) * h2(1);
4163 nrml[i * d1 + j * d2](1) =
4164 h1(2) * h2(0) - h1(0) * h2(2);
4165 nrml[i * d1 + j * d2](2) =
4166 h1(0) * h2(1) - h1(1) * h2(0);
4167
4168 // normalize Vector
4169 double norm =
4170 std::sqrt(std::pow(nrml[i * d1 + j * d2](0), 2.) +
4171 std::pow(nrml[i * d1 + j * d2](1), 2.) +
4172 std::pow(nrml[i * d1 + j * d2](2), 2.));
4173
4174 if (nrml[i * d1 + j * d2](1) < 0)
4175 norm *= -1.;
4176
4177 for (unsigned int k = 0; k < 3; ++k)
4178 nrml[i * d1 + j * d2](k) /= norm;
4179 }
4180 }
4181
4182 // setting up triangles
4183 for (unsigned int i = 0; i < n_subdivisions; ++i)
4184 for (unsigned int j = 0; j < n_subdivisions; ++j)
4185 {
4186 // down/left vertex of triangle
4187 const int dl = i * d1 + j * d2;
4188 if (flags.smooth)
4189 {
4190 // writing smooth_triangles
4191
4192 // down/right triangle
4193 out << "smooth_triangle {" << '\n'
4194 << "\t<" << ver[dl](0) << "," << patch.data(0, dl)
4195 << "," << ver[dl](1) << ">, <" << nrml[dl](0)
4196 << ", " << nrml[dl](1) << ", " << nrml[dl](2)
4197 << ">," << '\n';
4198 out << " \t<" << ver[dl + d1](0) << ","
4199 << patch.data(0, dl + d1) << "," << ver[dl + d1](1)
4200 << ">, <" << nrml[dl + d1](0) << ", "
4201 << nrml[dl + d1](1) << ", " << nrml[dl + d1](2)
4202 << ">," << '\n';
4203 out << "\t<" << ver[dl + d1 + d2](0) << ","
4204 << patch.data(0, dl + d1 + d2) << ","
4205 << ver[dl + d1 + d2](1) << ">, <"
4206 << nrml[dl + d1 + d2](0) << ", "
4207 << nrml[dl + d1 + d2](1) << ", "
4208 << nrml[dl + d1 + d2](2) << ">}" << '\n';
4209
4210 // upper/left triangle
4211 out << "smooth_triangle {" << '\n'
4212 << "\t<" << ver[dl](0) << "," << patch.data(0, dl)
4213 << "," << ver[dl](1) << ">, <" << nrml[dl](0)
4214 << ", " << nrml[dl](1) << ", " << nrml[dl](2)
4215 << ">," << '\n';
4216 out << "\t<" << ver[dl + d1 + d2](0) << ","
4217 << patch.data(0, dl + d1 + d2) << ","
4218 << ver[dl + d1 + d2](1) << ">, <"
4219 << nrml[dl + d1 + d2](0) << ", "
4220 << nrml[dl + d1 + d2](1) << ", "
4221 << nrml[dl + d1 + d2](2) << ">," << '\n';
4222 out << "\t<" << ver[dl + d2](0) << ","
4223 << patch.data(0, dl + d2) << "," << ver[dl + d2](1)
4224 << ">, <" << nrml[dl + d2](0) << ", "
4225 << nrml[dl + d2](1) << ", " << nrml[dl + d2](2)
4226 << ">}" << '\n';
4227 }
4228 else
4229 {
4230 // writing standard triangles down/right triangle
4231 out << "triangle {" << '\n'
4232 << "\t<" << ver[dl](0) << "," << patch.data(0, dl)
4233 << "," << ver[dl](1) << ">," << '\n';
4234 out << "\t<" << ver[dl + d1](0) << ","
4235 << patch.data(0, dl + d1) << "," << ver[dl + d1](1)
4236 << ">," << '\n';
4237 out << "\t<" << ver[dl + d1 + d2](0) << ","
4238 << patch.data(0, dl + d1 + d2) << ","
4239 << ver[dl + d1 + d2](1) << ">}" << '\n';
4240
4241 // upper/left triangle
4242 out << "triangle {" << '\n'
4243 << "\t<" << ver[dl](0) << "," << patch.data(0, dl)
4244 << "," << ver[dl](1) << ">," << '\n';
4245 out << "\t<" << ver[dl + d1 + d2](0) << ","
4246 << patch.data(0, dl + d1 + d2) << ","
4247 << ver[dl + d1 + d2](1) << ">," << '\n';
4248 out << "\t<" << ver[dl + d2](0) << ","
4249 << patch.data(0, dl + d2) << "," << ver[dl + d2](1)
4250 << ">}" << '\n';
4251 }
4252 }
4253 }
4254 else
4255 {
4256 // writing bicubic_patch
4257 Assert(n_subdivisions == 3,
4258 ExcDimensionMismatch(n_subdivisions, 3));
4259 out << '\n'
4260 << "bicubic_patch {" << '\n'
4261 << " type 0" << '\n'
4262 << " flatness 0" << '\n'
4263 << " u_steps 0" << '\n'
4264 << " v_steps 0" << '\n';
4265 for (int i = 0; i < 16; ++i)
4266 {
4267 out << "\t<" << ver[i](0) << "," << patch.data(0, i) << ","
4268 << ver[i](1) << ">";
4269 if (i != 15)
4270 out << ",";
4271 out << '\n';
4272 }
4273 out << " texture {Tex}" << '\n' << "}" << '\n';
4274 }
4275 }
4276
4277 if (!flags.bicubic_patch)
4278 {
4279 // the end of the mesh
4280 out << " texture {Tex}" << '\n' << "}" << '\n' << '\n';
4281 }
4282
4283 // make sure everything now gets to disk
4284 out.flush();
4285
4286 AssertThrow(out.fail() == false, ExcIO());
4287 }
4288 } // namespace
4289
4290
4291
4292 template <int dim, int spacedim>
4293 void
4295 const std::vector<Patch<dim, spacedim>> &patches,
4296 const std::vector<std::string> & data_names,
4297 const std::vector<
4298 std::tuple<unsigned int,
4299 unsigned int,
4300 std::string,
4302 const PovrayFlags &flags,
4303 std::ostream & out)
4304 {
4305 do_write_povray(patches, data_names, flags, out);
4306 }
4307
4308
4309
4310 template <int dim, int spacedim>
4311 void
4313 const std::vector<Patch<dim, spacedim>> & /*patches*/,
4314 const std::vector<std::string> & /*data_names*/,
4315 const std::vector<
4316 std::tuple<unsigned int,
4317 unsigned int,
4318 std::string,
4320 const EpsFlags & /*flags*/,
4321 std::ostream & /*out*/)
4322 {
4323 // not implemented, see the documentation of the function
4324 AssertThrow(dim == 2, ExcNotImplemented());
4325 }
4326
4327
4328 template <int spacedim>
4329 void
4331 const std::vector<Patch<2, spacedim>> &patches,
4332 const std::vector<std::string> & /*data_names*/,
4333 const std::vector<
4334 std::tuple<unsigned int,
4335 unsigned int,
4336 std::string,
4338 const EpsFlags &flags,
4339 std::ostream & out)
4340 {
4341 AssertThrow(out.fail() == false, ExcIO());
4342
4343#ifndef DEAL_II_WITH_MPI
4344 // verify that there are indeed patches to be written out. most of the
4345 // times, people just forget to call build_patches when there are no
4346 // patches, so a warning is in order. that said, the assertion is disabled
4347 // if we support MPI since then it can happen that on the coarsest mesh, a
4348 // processor simply has no cells it actually owns, and in that case it is
4349 // legit if there are no patches
4350 Assert(patches.size() > 0, ExcNoPatches());
4351#else
4352 if (patches.size() == 0)
4353 return;
4354#endif
4355
4356 // set up an array of cells to be written later. this array holds the cells
4357 // of all the patches as projected to the plane perpendicular to the line of
4358 // sight.
4359 //
4360 // note that they are kept sorted by the set, where we chose the value of
4361 // the center point of the cell along the line of sight as value for sorting
4362 std::multiset<EpsCell2d> cells;
4363
4364 // two variables in which we will store the minimum and maximum values of
4365 // the field to be used for colorization
4366 float min_color_value = std::numeric_limits<float>::max();
4367 float max_color_value = std::numeric_limits<float>::min();
4368
4369 // Array for z-coordinates of points. The elevation determined by a function
4370 // if spacedim=2 or the z-coordinate of the grid point if spacedim=3
4371 double heights[4] = {0, 0, 0, 0};
4372
4373 // compute the cells for output and enter them into the set above note that
4374 // since dim==2, we have exactly four vertices per patch and per cell
4375 for (const auto &patch : patches)
4376 {
4377 const unsigned int n_subdivisions = patch.n_subdivisions;
4378 const unsigned int n = n_subdivisions + 1;
4379 const unsigned int d1 = 1;
4380 const unsigned int d2 = n;
4381
4382 for (unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
4383 for (unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
4384 {
4385 Point<spacedim> points[4];
4386 points[0] =
4387 get_equispaced_location(patch, {i1, i2}, n_subdivisions);
4388 points[1] =
4389 get_equispaced_location(patch, {i1 + 1, i2}, n_subdivisions);
4390 points[2] =
4391 get_equispaced_location(patch, {i1, i2 + 1}, n_subdivisions);
4392 points[3] = get_equispaced_location(patch,
4393 {i1 + 1, i2 + 1},
4394 n_subdivisions);
4395
4396 switch (spacedim)
4397 {
4398 case 2:
4399 Assert((flags.height_vector < patch.data.n_rows()) ||
4400 patch.data.n_rows() == 0,
4402 0,
4403 patch.data.n_rows()));
4404 heights[0] =
4405 patch.data.n_rows() != 0 ?
4406 patch.data(flags.height_vector, i1 * d1 + i2 * d2) *
4407 flags.z_scaling :
4408 0;
4409 heights[1] = patch.data.n_rows() != 0 ?
4410 patch.data(flags.height_vector,
4411 (i1 + 1) * d1 + i2 * d2) *
4412 flags.z_scaling :
4413 0;
4414 heights[2] = patch.data.n_rows() != 0 ?
4415 patch.data(flags.height_vector,
4416 i1 * d1 + (i2 + 1) * d2) *
4417 flags.z_scaling :
4418 0;
4419 heights[3] = patch.data.n_rows() != 0 ?
4420 patch.data(flags.height_vector,
4421 (i1 + 1) * d1 + (i2 + 1) * d2) *
4422 flags.z_scaling :
4423 0;
4424
4425 break;
4426 case 3:
4427 // Copy z-coordinates into the height vector
4428 for (unsigned int i = 0; i < 4; ++i)
4429 heights[i] = points[i](2);
4430 break;
4431 default:
4432 Assert(false, ExcNotImplemented());
4433 }
4434
4435
4436 // now compute the projection of the bilinear cell given by the
4437 // four vertices and their heights and write them to a proper cell
4438 // object. note that we only need the first two components of the
4439 // projected position for output, but we need the value along the
4440 // line of sight for sorting the cells for back-to- front-output
4441 //
4442 // this computation was first written by Stefan Nauber. please
4443 // no-one ask me why it works that way (or may be not), especially
4444 // not about the angles and the sign of the height field, I don't
4445 // know it.
4446 EpsCell2d eps_cell;
4447 const double pi = numbers::PI;
4448 const double cx =
4449 -std::cos(pi - flags.azimut_angle * 2 * pi / 360.),
4450 cz = -std::cos(flags.turn_angle * 2 * pi / 360.),
4451 sx =
4452 std::sin(pi - flags.azimut_angle * 2 * pi / 360.),
4453 sz = std::sin(flags.turn_angle * 2 * pi / 360.);
4454 for (unsigned int vertex = 0; vertex < 4; ++vertex)
4455 {
4456 const double x = points[vertex](0), y = points[vertex](1),
4457 z = -heights[vertex];
4458
4459 eps_cell.vertices[vertex](0) = -cz * x + sz * y;
4460 eps_cell.vertices[vertex](1) =
4461 -cx * sz * x - cx * cz * y - sx * z;
4462
4463 // ( 1 0 0 )
4464 // D1 = ( 0 cx -sx )
4465 // ( 0 sx cx )
4466
4467 // ( cy 0 sy )
4468 // Dy = ( 0 1 0 )
4469 // (-sy 0 cy )
4470
4471 // ( cz -sz 0 )
4472 // Dz = ( sz cz 0 )
4473 // ( 0 0 1 )
4474
4475 // ( cz -sz 0 )( 1 0 0 )(x) (
4476 // cz*x-sz*(cx*y-sx*z)+0*(sx*y+cx*z) )
4477 // Dxz = ( sz cz 0 )( 0 cx -sx )(y) = (
4478 // sz*x+cz*(cx*y-sx*z)+0*(sx*y+cx*z) )
4479 // ( 0 0 1 )( 0 sx cx )(z) ( 0*x+
4480 // *(cx*y-sx*z)+1*(sx*y+cx*z) )
4481 }
4482
4483 // compute coordinates of center of cell
4484 const Point<spacedim> center_point =
4485 (points[0] + points[1] + points[2] + points[3]) / 4;
4486 const double center_height =
4487 -(heights[0] + heights[1] + heights[2] + heights[3]) / 4;
4488
4489 // compute the depth into the picture
4490 eps_cell.depth = -sx * sz * center_point(0) -
4491 sx * cz * center_point(1) + cx * center_height;
4492
4493 if (flags.draw_cells && flags.shade_cells)
4494 {
4495 Assert((flags.color_vector < patch.data.n_rows()) ||
4496 patch.data.n_rows() == 0,
4498 0,
4499 patch.data.n_rows()));
4500 const double color_values[4] = {
4501 patch.data.n_rows() != 0 ?
4502 patch.data(flags.color_vector, i1 * d1 + i2 * d2) :
4503 1,
4504
4505 patch.data.n_rows() != 0 ?
4506 patch.data(flags.color_vector, (i1 + 1) * d1 + i2 * d2) :
4507 1,
4508
4509 patch.data.n_rows() != 0 ?
4510 patch.data(flags.color_vector, i1 * d1 + (i2 + 1) * d2) :
4511 1,
4512
4513 patch.data.n_rows() != 0 ?
4514 patch.data(flags.color_vector,
4515 (i1 + 1) * d1 + (i2 + 1) * d2) :
4516 1};
4517
4518 // set color value to average of the value at the vertices
4519 eps_cell.color_value = (color_values[0] + color_values[1] +
4520 color_values[3] + color_values[2]) /
4521 4;
4522
4523 // update bounds of color field
4524 min_color_value =
4525 std::min(min_color_value, eps_cell.color_value);
4526 max_color_value =
4527 std::max(max_color_value, eps_cell.color_value);
4528 }
4529
4530 // finally add this cell
4531 cells.insert(eps_cell);
4532 }
4533 }
4534
4535 // find out minimum and maximum x and y coordinates to compute offsets and
4536 // scaling factors
4537 double x_min = cells.begin()->vertices[0](0);
4538 double x_max = x_min;
4539 double y_min = cells.begin()->vertices[0](1);
4540 double y_max = y_min;
4541
4542 for (const auto &cell : cells)
4543 for (const auto &vertex : cell.vertices)
4544 {
4545 x_min = std::min(x_min, vertex(0));
4546 x_max = std::max(x_max, vertex(0));
4547 y_min = std::min(y_min, vertex(1));
4548 y_max = std::max(y_max, vertex(1));
4549 }
4550
4551 // scale in x-direction such that in the output 0 <= x <= 300. don't scale
4552 // in y-direction to preserve the shape of the triangulation
4553 const double scale =
4554 (flags.size /
4555 (flags.size_type == EpsFlags::width ? x_max - x_min : y_min - y_max));
4556
4557 const Point<2> offset(x_min, y_min);
4558
4559
4560 // now write preamble
4561 {
4562 out << "%!PS-Adobe-2.0 EPSF-1.2" << '\n'
4563 << "%%Title: deal.II Output" << '\n'
4564 << "%%Creator: the deal.II library" << '\n'
4565 << "%%Creation Date: " << Utilities::System::get_date() << " - "
4566 << Utilities::System::get_time() << '\n'
4567 << "%%BoundingBox: "
4568 // lower left corner
4569 << "0 0 "
4570 // upper right corner
4571 << static_cast<unsigned int>((x_max - x_min) * scale + 0.5) << ' '
4572 << static_cast<unsigned int>((y_max - y_min) * scale + 0.5) << '\n';
4573
4574 // define some abbreviations to keep the output small:
4575 // m=move turtle to
4576 // l=define a line
4577 // s=set rgb color
4578 // sg=set gray value
4579 // lx=close the line and plot the line
4580 // lf=close the line and fill the interior
4581 out << "/m {moveto} bind def" << '\n'
4582 << "/l {lineto} bind def" << '\n'
4583 << "/s {setrgbcolor} bind def" << '\n'
4584 << "/sg {setgray} bind def" << '\n'
4585 << "/lx {lineto closepath stroke} bind def" << '\n'
4586 << "/lf {lineto closepath fill} bind def" << '\n';
4587
4588 out << "%%EndProlog" << '\n' << '\n';
4589 // set fine lines
4590 out << flags.line_width << " setlinewidth" << '\n';
4591 }
4592
4593 // check if min and max values for the color are actually different. If
4594 // that is not the case (such things happen, for example, in the very first
4595 // time step of a time dependent problem, if the initial values are zero),
4596 // all values are equal, and then we can draw everything in an arbitrary
4597 // color. Thus, change one of the two values arbitrarily
4598 if (max_color_value == min_color_value)
4599 max_color_value = min_color_value + 1;
4600
4601 // now we've got all the information we need. write the cells. note: due to
4602 // the ordering, we traverse the list of cells back-to-front
4603 for (const auto &cell : cells)
4604 {
4605 if (flags.draw_cells)
4606 {
4607 if (flags.shade_cells)
4608 {
4609 const EpsFlags::RgbValues rgb_values =
4610 (*flags.color_function)(cell.color_value,
4611 min_color_value,
4612 max_color_value);
4613
4614 // write out color
4615 if (rgb_values.is_grey())
4616 out << rgb_values.red << " sg ";
4617 else
4618 out << rgb_values.red << ' ' << rgb_values.green << ' '
4619 << rgb_values.blue << " s ";
4620 }
4621 else
4622 out << "1 sg ";
4623
4624 out << (cell.vertices[0] - offset) * scale << " m "
4625 << (cell.vertices[1] - offset) * scale << " l "
4626 << (cell.vertices[3] - offset) * scale << " l "
4627 << (cell.vertices[2] - offset) * scale << " lf" << '\n';
4628 }
4629
4630 if (flags.draw_mesh)
4631 out << "0 sg " // draw lines in black
4632 << (cell.vertices[0] - offset) * scale << " m "
4633 << (cell.vertices[1] - offset) * scale << " l "
4634 << (cell.vertices[3] - offset) * scale << " l "
4635 << (cell.vertices[2] - offset) * scale << " lx" << '\n';
4636 }
4637 out << "showpage" << '\n';
4638
4639 out.flush();
4640
4641 AssertThrow(out.fail() == false, ExcIO());
4642 }
4643
4644
4645
4646 template <int dim, int spacedim>
4647 void
4649 const std::vector<Patch<dim, spacedim>> &patches,
4650 const std::vector<std::string> & data_names,
4651 const std::vector<
4652 std::tuple<unsigned int,
4653 unsigned int,
4654 std::string,
4656 const GmvFlags &flags,
4657 std::ostream & out)
4658 {
4659 // The gmv format does not support cells that only consist of a single
4660 // point. It does support the output of point data using the keyword
4661 // 'tracers' instead of 'nodes' and 'cells', but this output format is
4662 // currently not implemented.
4663 AssertThrow(dim > 0, ExcNotImplemented());
4664
4665 Assert(dim <= 3, ExcNotImplemented());
4666 AssertThrow(out.fail() == false, ExcIO());
4667
4668#ifndef DEAL_II_WITH_MPI
4669 // verify that there are indeed patches to be written out. most of the
4670 // times, people just forget to call build_patches when there are no
4671 // patches, so a warning is in order. that said, the assertion is disabled
4672 // if we support MPI since then it can happen that on the coarsest mesh, a
4673 // processor simply has no cells it actually owns, and in that case it is
4674 // legit if there are no patches
4675 Assert(patches.size() > 0, ExcNoPatches());
4676#else
4677 if (patches.size() == 0)
4678 return;
4679#endif
4680
4681 GmvStream gmv_out(out, flags);
4682 const unsigned int n_data_sets = data_names.size();
4683 // check against # of data sets in first patch. checks against all other
4684 // patches are made in write_gmv_reorder_data_vectors
4685 Assert((patches[0].data.n_rows() == n_data_sets &&
4686 !patches[0].points_are_available) ||
4687 (patches[0].data.n_rows() == n_data_sets + spacedim &&
4688 patches[0].points_are_available),
4689 ExcDimensionMismatch(patches[0].points_are_available ?
4690 (n_data_sets + spacedim) :
4691 n_data_sets,
4692 patches[0].data.n_rows()));
4693
4694 //---------------------
4695 // preamble
4696 out << "gmvinput ascii" << '\n' << '\n';
4697
4698 // first count the number of cells and cells for later use
4699 unsigned int n_nodes;
4700 unsigned int n_cells;
4701 std::tie(n_nodes, n_cells) = count_nodes_and_cells(patches);
4702
4703 // For the format we write here, we need to write all node values relating
4704 // to one variable at a time. We could in principle do this by looping
4705 // over all patches and extracting the values corresponding to the one
4706 // variable we're dealing with right now, and then start the process over
4707 // for the next variable with another loop over all patches.
4708 //
4709 // An easier way is to create a global table that for each variable
4710 // lists all values. This copying of data vectors can be done in the
4711 // background while we're already working on vertices and cells,
4712 // so do this on a separate task and when wanting to write out the
4713 // data, we wait for that task to finish.
4715 create_global_data_table_task = Threads::new_task(
4716 [&patches]() { return create_global_data_table(patches); });
4717
4718 //-----------------------------
4719 // first make up a list of used vertices along with their coordinates
4720 //
4721 // note that we have to print 3 dimensions
4722 out << "nodes " << n_nodes << '\n';
4723 for (unsigned int d = 0; d < spacedim; ++d)
4724 {
4725 gmv_out.selected_component = d;
4726 write_nodes(patches, gmv_out);
4727 out << '\n';
4728 }
4729 gmv_out.selected_component = numbers::invalid_unsigned_int;
4730
4731 for (unsigned int d = spacedim; d < 3; ++d)
4732 {
4733 for (unsigned int i = 0; i < n_nodes; ++i)
4734 out << "0 ";
4735 out << '\n';
4736 }
4737
4738 //-------------------------------
4739 // now for the cells. note that vertices are counted from 1 onwards
4740 out << "cells " << n_cells << '\n';
4741 write_cells(patches, gmv_out);
4742
4743 //-------------------------------------
4744 // data output.
4745 out << "variable" << '\n';
4746
4747 // Wait for the reordering to be done and retrieve the reordered data:
4748 const Table<2, double> data_vectors =
4749 std::move(*create_global_data_table_task.return_value());
4750
4751 // then write data. the '1' means: node data (as opposed to cell data, which
4752 // we do not support explicitly here)
4753 for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
4754 {
4755 out << data_names[data_set] << " 1" << '\n';
4756 std::copy(data_vectors[data_set].begin(),
4757 data_vectors[data_set].end(),
4758 std::ostream_iterator<double>(out, " "));
4759 out << '\n' << '\n';
4760 }
4761
4762
4763
4764 // end of variable section
4765 out << "endvars" << '\n';
4766
4767 // end of output
4768 out << "endgmv" << '\n';
4769
4770 // make sure everything now gets to disk
4771 out.flush();
4772
4773 // assert the stream is still ok
4774 AssertThrow(out.fail() == false, ExcIO());
4775 }
4776
4777
4778
4779 template <int dim, int spacedim>
4780 void
4782 const std::vector<Patch<dim, spacedim>> &patches,
4783 const std::vector<std::string> & data_names,
4784 const std::vector<
4785 std::tuple<unsigned int,
4786 unsigned int,
4787 std::string,
4789 const TecplotFlags &flags,
4790 std::ostream & out)
4791 {
4792 AssertThrow(out.fail() == false, ExcIO());
4793
4794 // The FEBLOCK or FEPOINT formats of tecplot only allows full elements (e.g.
4795 // triangles), not single points. Other tecplot format allow point output,
4796 // but they are currently not implemented.
4797 AssertThrow(dim > 0, ExcNotImplemented());
4798
4799#ifndef DEAL_II_WITH_MPI
4800 // verify that there are indeed patches to be written out. most of the
4801 // times, people just forget to call build_patches when there are no
4802 // patches, so a warning is in order. that said, the assertion is disabled
4803 // if we support MPI since then it can happen that on the coarsest mesh, a
4804 // processor simply has no cells it actually owns, and in that case it is
4805 // legit if there are no patches
4806 Assert(patches.size() > 0, ExcNoPatches());
4807#else
4808 if (patches.size() == 0)
4809 return;
4810#endif
4811
4812 TecplotStream tecplot_out(out, flags);
4813
4814 const unsigned int n_data_sets = data_names.size();
4815 // check against # of data sets in first patch. checks against all other
4816 // patches are made in write_gmv_reorder_data_vectors
4817 Assert((patches[0].data.n_rows() == n_data_sets &&
4818 !patches[0].points_are_available) ||
4819 (patches[0].data.n_rows() == n_data_sets + spacedim &&
4820 patches[0].points_are_available),
4821 ExcDimensionMismatch(patches[0].points_are_available ?
4822 (n_data_sets + spacedim) :
4823 n_data_sets,
4824 patches[0].data.n_rows()));
4825
4826 // first count the number of cells and cells for later use
4827 unsigned int n_nodes;
4828 unsigned int n_cells;
4829 std::tie(n_nodes, n_cells) = count_nodes_and_cells(patches);
4830
4831 //---------
4832 // preamble
4833 {
4834 out
4835 << "# This file was generated by the deal.II library." << '\n'
4836 << "# Date = " << Utilities::System::get_date() << '\n'
4837 << "# Time = " << Utilities::System::get_time() << '\n'
4838 << "#" << '\n'
4839 << "# For a description of the Tecplot format see the Tecplot documentation."
4840 << '\n'
4841 << "#" << '\n';
4842
4843
4844 out << "Variables=";
4845
4846 switch (spacedim)
4847 {
4848 case 1:
4849 out << "\"x\"";
4850 break;
4851 case 2:
4852 out << "\"x\", \"y\"";
4853 break;
4854 case 3:
4855 out << "\"x\", \"y\", \"z\"";
4856 break;
4857 default:
4858 Assert(false, ExcNotImplemented());
4859 }
4860
4861 for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
4862 out << ", \"" << data_names[data_set] << "\"";
4863
4864 out << '\n';
4865
4866 out << "zone ";
4867 if (flags.zone_name)
4868 out << "t=\"" << flags.zone_name << "\" ";
4869
4870 if (flags.solution_time >= 0.0)
4871 out << "strandid=1, solutiontime=" << flags.solution_time << ", ";
4872
4873 out << "f=feblock, n=" << n_nodes << ", e=" << n_cells
4874 << ", et=" << tecplot_cell_type[dim] << '\n';
4875 }
4876
4877
4878 // For the format we write here, we need to write all node values relating
4879 // to one variable at a time. We could in principle do this by looping
4880 // over all patches and extracting the values corresponding to the one
4881 // variable we're dealing with right now, and then start the process over
4882 // for the next variable with another loop over all patches.
4883 //
4884 // An easier way is to create a global table that for each variable
4885 // lists all values. This copying of data vectors can be done in the
4886 // background while we're already working on vertices and cells,
4887 // so do this on a separate task and when wanting to write out the
4888 // data, we wait for that task to finish.
4890 create_global_data_table_task = Threads::new_task(
4891 [&patches]() { return create_global_data_table(patches); });
4892
4893 //-----------------------------
4894 // first make up a list of used vertices along with their coordinates
4895
4896
4897 for (unsigned int d = 0; d < spacedim; ++d)
4898 {
4899 tecplot_out.selected_component = d;
4900 write_nodes(patches, tecplot_out);
4901 out << '\n';
4902 }
4903
4904
4905 //-------------------------------------
4906 // data output.
4907 //
4908 // Wait for the reordering to be done and retrieve the reordered data:
4909 const Table<2, double> data_vectors =
4910 std::move(*create_global_data_table_task.return_value());
4911
4912 // then write data.
4913 for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
4914 {
4915 std::copy(data_vectors[data_set].begin(),
4916 data_vectors[data_set].end(),
4917 std::ostream_iterator<double>(out, "\n"));
4918 out << '\n';
4919 }
4920
4921 write_cells(patches, tecplot_out);
4922
4923 // make sure everything now gets to disk
4924 out.flush();
4925
4926 // assert the stream is still ok
4927 AssertThrow(out.fail() == false, ExcIO());
4928 }
4929
4930
4931
4932 template <int dim, int spacedim>
4933 void
4935 const std::vector<Patch<dim, spacedim>> &patches,
4936 const std::vector<std::string> & data_names,
4937 const std::vector<
4938 std::tuple<unsigned int,
4939 unsigned int,
4940 std::string,
4942 & nonscalar_data_ranges,
4943 const VtkFlags &flags,
4944 std::ostream & out)
4945 {
4946 AssertThrow(out.fail() == false, ExcIO());
4947
4948#ifndef DEAL_II_WITH_MPI
4949 // verify that there are indeed patches to be written out. most of the
4950 // times, people just forget to call build_patches when there are no
4951 // patches, so a warning is in order. that said, the assertion is disabled
4952 // if we support MPI since then it can happen that on the coarsest mesh, a
4953 // processor simply has no cells it actually owns, and in that case it is
4954 // legit if there are no patches
4955 Assert(patches.size() > 0, ExcNoPatches());
4956#else
4957 if (patches.size() == 0)
4958 return;
4959#endif
4960
4961 VtkStream vtk_out(out, flags);
4962
4963 const unsigned int n_data_sets = data_names.size();
4964 // check against # of data sets in first patch.
4965 if (patches[0].points_are_available)
4966 {
4967 AssertDimension(n_data_sets + spacedim, patches[0].data.n_rows())
4968 }
4969 else
4970 {
4971 AssertDimension(n_data_sets, patches[0].data.n_rows())
4972 }
4973
4974 //---------------------
4975 // preamble
4976 {
4977 out << "# vtk DataFile Version 3.0" << '\n'
4978 << "#This file was generated by the deal.II library";
4979 if (flags.print_date_and_time)
4980 {
4981 out << " on " << Utilities::System::get_date() << " at "
4983 }
4984 else
4985 out << '.';
4986 out << '\n' << "ASCII" << '\n';
4987 // now output the data header
4988 out << "DATASET UNSTRUCTURED_GRID\n" << '\n';
4989 }
4990
4991 // if desired, output time and cycle of the simulation, following the
4992 // instructions at
4993 // http://www.visitusers.org/index.php?title=Time_and_Cycle_in_VTK_files
4994 {
4995 const unsigned int n_metadata =
4996 ((flags.cycle != std::numeric_limits<unsigned int>::min() ? 1 : 0) +
4997 (flags.time != std::numeric_limits<double>::min() ? 1 : 0));
4998 if (n_metadata > 0)
4999 {
5000 out << "FIELD FieldData " << n_metadata << '\n';
5001
5002 if (flags.cycle != std::numeric_limits<unsigned int>::min())
5003 {
5004 out << "CYCLE 1 1 int\n" << flags.cycle << '\n';
5005 }
5006 if (flags.time != std::numeric_limits<double>::min())
5007 {
5008 out << "TIME 1 1 double\n" << flags.time << '\n';
5009 }
5010 }
5011 }
5012
5013 // first count the number of cells and cells for later use
5014 unsigned int n_nodes;
5015 unsigned int n_cells;
5016 unsigned int n_points_and_n_cells;
5017 std::tie(n_nodes, n_cells, n_points_and_n_cells) =
5018 count_nodes_and_cells_and_points(patches, flags.write_higher_order_cells);
5019
5020 // For the format we write here, we need to write all node values relating
5021 // to one variable at a time. We could in principle do this by looping
5022 // over all patches and extracting the values corresponding to the one
5023 // variable we're dealing with right now, and then start the process over
5024 // for the next variable with another loop over all patches.
5025 //
5026 // An easier way is to create a global table that for each variable
5027 // lists all values. This copying of data vectors can be done in the
5028 // background while we're already working on vertices and cells,
5029 // so do this on a separate task and when wanting to write out the
5030 // data, we wait for that task to finish.
5032 create_global_data_table_task = Threads::new_task(
5033 [&patches]() { return create_global_data_table(patches); });
5034
5035 //-----------------------------
5036 // first make up a list of used vertices along with their coordinates
5037 //
5038 // note that we have to print d=1..3 dimensions
5039 out << "POINTS " << n_nodes << " double" << '\n';
5040 write_nodes(patches, vtk_out);
5041 out << '\n';
5042 //-------------------------------
5043 // now for the cells
5044 out << "CELLS " << n_cells << ' ' << n_points_and_n_cells << '\n';
5045 if (flags.write_higher_order_cells)
5046 write_high_order_cells(patches, vtk_out, /* legacy_format = */ true);
5047 else
5048 write_cells(patches, vtk_out);
5049 out << '\n';
5050 // next output the types of the cells. since all cells are the same, this is
5051 // simple
5052 out << "CELL_TYPES " << n_cells << '\n';
5053
5054 // need to distinguish between linear cells, simplex cells (linear or
5055 // quadratic), and high order cells
5056 for (const auto &patch : patches)
5057 {
5058 const auto vtk_cell_id =
5059 extract_vtk_patch_info(patch, flags.write_higher_order_cells);
5060
5061 for (unsigned int i = 0; i < vtk_cell_id[1]; ++i)
5062 out << ' ' << vtk_cell_id[0];
5063 }
5064
5065 out << '\n';
5066 //-------------------------------------
5067 // data output.
5068
5069 // Wait for the reordering to be done and retrieve the reordered data:
5070 const Table<2, double> data_vectors =
5071 std::move(*create_global_data_table_task.return_value());
5072
5073 // then write data. the 'POINT_DATA' means: node data (as opposed to cell
5074 // data, which we do not support explicitly here). all following data sets
5075 // are point data
5076 out << "POINT_DATA " << n_nodes << '\n';
5077
5078 // when writing, first write out all vector data, then handle the scalar
5079 // data sets that have been left over
5080 std::vector<bool> data_set_written(n_data_sets, false);
5081 for (const auto &nonscalar_data_range : nonscalar_data_ranges)
5082 {
5083 AssertThrow(std::get<3>(nonscalar_data_range) !=
5085 ExcMessage(
5086 "The VTK writer does not currently support outputting "
5087 "tensor data. Use the VTU writer instead."));
5088
5089 AssertThrow(std::get<1>(nonscalar_data_range) >=
5090 std::get<0>(nonscalar_data_range),
5091 ExcLowerRange(std::get<1>(nonscalar_data_range),
5092 std::get<0>(nonscalar_data_range)));
5093 AssertThrow(std::get<1>(nonscalar_data_range) < n_data_sets,
5094 ExcIndexRange(std::get<1>(nonscalar_data_range),
5095 0,
5096 n_data_sets));
5097 AssertThrow(std::get<1>(nonscalar_data_range) + 1 -
5098 std::get<0>(nonscalar_data_range) <=
5099 3,
5100 ExcMessage(
5101 "Can't declare a vector with more than 3 components "
5102 "in VTK"));
5103
5104 // mark these components as already written:
5105 for (unsigned int i = std::get<0>(nonscalar_data_range);
5106 i <= std::get<1>(nonscalar_data_range);
5107 ++i)
5108 data_set_written[i] = true;
5109
5110 // write the header. concatenate all the component names with double
5111 // underscores unless a vector name has been specified
5112 out << "VECTORS ";
5113
5114 if (!std::get<2>(nonscalar_data_range).empty())
5115 out << std::get<2>(nonscalar_data_range);
5116 else
5117 {
5118 for (unsigned int i = std::get<0>(nonscalar_data_range);
5119 i < std::get<1>(nonscalar_data_range);
5120 ++i)
5121 out << data_names[i] << "__";
5122 out << data_names[std::get<1>(nonscalar_data_range)];
5123 }
5124
5125 out << " double" << '\n';
5126
5127 // now write data. pad all vectors to have three components
5128 for (unsigned int n = 0; n < n_nodes; ++n)
5129 {
5130 switch (std::get<1>(nonscalar_data_range) -
5131 std::get<0>(nonscalar_data_range))
5132 {
5133 case 0:
5134 out << data_vectors(std::get<0>(nonscalar_data_range), n)
5135 << " 0 0" << '\n';
5136 break;
5137
5138 case 1:
5139 out << data_vectors(std::get<0>(nonscalar_data_range), n)
5140 << ' '
5141 << data_vectors(std::get<0>(nonscalar_data_range) + 1, n)
5142 << " 0" << '\n';
5143 break;
5144 case 2:
5145 out << data_vectors(std::get<0>(nonscalar_data_range), n)
5146 << ' '
5147 << data_vectors(std::get<0>(nonscalar_data_range) + 1, n)
5148 << ' '
5149 << data_vectors(std::get<0>(nonscalar_data_range) + 2, n)
5150 << '\n';
5151 break;
5152
5153 default:
5154 // VTK doesn't support anything else than vectors with 1, 2,
5155 // or 3 components
5156 Assert(false, ExcInternalError());
5157 }
5158 }
5159 }
5160
5161 // now do the left over scalar data sets
5162 for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
5163 if (data_set_written[data_set] == false)
5164 {
5165 out << "SCALARS " << data_names[data_set] << " double 1" << '\n'
5166 << "LOOKUP_TABLE default" << '\n';
5167 std::copy(data_vectors[data_set].begin(),
5168 data_vectors[data_set].end(),
5169 std::ostream_iterator<double>(out, " "));
5170 out << '\n';
5171 }
5172
5173 // make sure everything now gets to disk
5174 out.flush();
5175
5176 // assert the stream is still ok
5177 AssertThrow(out.fail() == false, ExcIO());
5178 }
5179
5180
5181 void
5182 write_vtu_header(std::ostream &out, const VtkFlags &flags)
5183 {
5184 AssertThrow(out.fail() == false, ExcIO());
5185 out << "<?xml version=\"1.0\" ?> \n";
5186 out << "<!-- \n";
5187 out << "# vtk DataFile Version 3.0" << '\n'
5188 << "#This file was generated by the deal.II library";
5189 if (flags.print_date_and_time)
5190 {
5191 out << " on " << Utilities::System::get_time() << " at "
5193 }
5194 else
5195 out << '.';
5196 out << "\n-->\n";
5197
5198 if (flags.write_higher_order_cells)
5199 out << "<VTKFile type=\"UnstructuredGrid\" version=\"2.2\"";
5200 else
5201 out << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\"";
5202 if (deal_ii_with_zlib &&
5204 out << " compressor=\"vtkZLibDataCompressor\"";
5205#ifdef DEAL_II_WORDS_BIGENDIAN
5206 out << " byte_order=\"BigEndian\"";
5207#else
5208 out << " byte_order=\"LittleEndian\"";
5209#endif
5210 out << ">";
5211 out << '\n';
5212 out << "<UnstructuredGrid>";
5213 out << '\n';
5214 }
5215
5216
5217
5218 void
5219 write_vtu_footer(std::ostream &out)
5220 {
5221 AssertThrow(out.fail() == false, ExcIO());
5222 out << " </UnstructuredGrid>\n";
5223 out << "</VTKFile>\n";
5224 }
5225
5226
5227
5228 template <int dim, int spacedim>
5229 void
5231 const std::vector<Patch<dim, spacedim>> &patches,
5232 const std::vector<std::string> & data_names,
5233 const std::vector<
5234 std::tuple<unsigned int,
5235 unsigned int,
5236 std::string,
5238 & nonscalar_data_ranges,
5239 const VtkFlags &flags,
5240 std::ostream & out)
5241 {
5242 write_vtu_header(out, flags);
5243 write_vtu_main(patches, data_names, nonscalar_data_ranges, flags, out);
5244 write_vtu_footer(out);
5245
5246 out << std::flush;
5247 }
5248
5249
5250 template <int dim, int spacedim>
5251 void
5253 const std::vector<Patch<dim, spacedim>> &patches,
5254 const std::vector<std::string> & data_names,
5255 const std::vector<
5256 std::tuple<unsigned int,
5257 unsigned int,
5258 std::string,
5260 & nonscalar_data_ranges,
5261 const VtkFlags &flags,
5262 std::ostream & out)
5263 {
5264 AssertThrow(out.fail() == false, ExcIO());
5265
5266 // If the user provided physical units, make sure that they don't contain
5267 // quote characters as this would make the VTU file invalid XML and
5268 // probably lead to all sorts of difficult error messages. Other than that,
5269 // trust the user that whatever they provide makes sense somehow.
5270 for (const auto &unit : flags.physical_units)
5271 {
5272 (void)unit;
5273 Assert(
5274 unit.second.find('\"') == std::string::npos,
5275 ExcMessage(
5276 "A physical unit you provided, <" + unit.second +
5277 ">, contained a quotation mark character. This is not allowed."));
5278 }
5279
5280#ifndef DEAL_II_WITH_MPI
5281 // verify that there are indeed patches to be written out. most of the
5282 // times, people just forget to call build_patches when there are no
5283 // patches, so a warning is in order. that said, the assertion is disabled
5284 // if we support MPI since then it can happen that on the coarsest mesh, a
5285 // processor simply has no cells it actually owns, and in that case it is
5286 // legit if there are no patches
5287 Assert(patches.size() > 0, ExcNoPatches());
5288#else
5289 if (patches.size() == 0)
5290 {
5291 // we still need to output a valid vtu file, because other CPUs might
5292 // output data. This is the minimal file that is accepted by paraview
5293 // and visit. if we remove the field definitions, visit is complaining.
5294 out << "<Piece NumberOfPoints=\"0\" NumberOfCells=\"0\" >\n"
5295 << "<Cells>\n"
5296 << "<DataArray type=\"UInt8\" Name=\"types\"></DataArray>\n"
5297 << "</Cells>\n"
5298 << " <PointData Scalars=\"scalars\">\n";
5299 std::vector<bool> data_set_written(data_names.size(), false);
5300 for (const auto &nonscalar_data_range : nonscalar_data_ranges)
5301 {
5302 // mark these components as already written:
5303 for (unsigned int i = std::get<0>(nonscalar_data_range);
5304 i <= std::get<1>(nonscalar_data_range);
5305 ++i)
5306 data_set_written[i] = true;
5307
5308 // write the header. concatenate all the component names with double
5309 // underscores unless a vector name has been specified
5310 out << " <DataArray type=\"Float32\" Name=\"";
5311
5312 if (!std::get<2>(nonscalar_data_range).empty())
5313 out << std::get<2>(nonscalar_data_range);
5314 else
5315 {
5316 for (unsigned int i = std::get<0>(nonscalar_data_range);
5317 i < std::get<1>(nonscalar_data_range);
5318 ++i)
5319 out << data_names[i] << "__";
5320 out << data_names[std::get<1>(nonscalar_data_range)];
5321 }
5322
5323 out << "\" NumberOfComponents=\"3\"></DataArray>\n";
5324 }
5325
5326 for (unsigned int data_set = 0; data_set < data_names.size();
5327 ++data_set)
5328 if (data_set_written[data_set] == false)
5329 {
5330 out << " <DataArray type=\"Float32\" Name=\""
5331 << data_names[data_set] << "\"></DataArray>\n";
5332 }
5333
5334 out << " </PointData>\n";
5335 out << "</Piece>\n";
5336
5337 out << std::flush;
5338
5339 return;
5340 }
5341#endif
5342
5343 // first up: metadata
5344 //
5345 // if desired, output time and cycle of the simulation, following the
5346 // instructions at
5347 // http://www.visitusers.org/index.php?title=Time_and_Cycle_in_VTK_files
5348 {
5349 const unsigned int n_metadata =
5350 ((flags.cycle != std::numeric_limits<unsigned int>::min() ? 1 : 0) +
5351 (flags.time != std::numeric_limits<double>::min() ? 1 : 0));
5352 if (n_metadata > 0)
5353 out << "<FieldData>\n";
5354
5355 if (flags.cycle != std::numeric_limits<unsigned int>::min())
5356 {
5357 out
5358 << "<DataArray type=\"Float32\" Name=\"CYCLE\" NumberOfTuples=\"1\" format=\"ascii\">"
5359 << flags.cycle << "</DataArray>\n";
5360 }
5361 if (flags.time != std::numeric_limits<double>::min())
5362 {
5363 out
5364 << "<DataArray type=\"Float32\" Name=\"TIME\" NumberOfTuples=\"1\" format=\"ascii\">"
5365 << flags.time << "</DataArray>\n";
5366 }
5367
5368 if (n_metadata > 0)
5369 out << "</FieldData>\n";
5370 }
5371
5372
5373 const unsigned int n_data_sets = data_names.size();
5374 // check against # of data sets in first patch. checks against all other
5375 // patches are made in write_gmv_reorder_data_vectors
5376 if (patches[0].points_are_available)
5377 {
5378 AssertDimension(n_data_sets + spacedim, patches[0].data.n_rows())
5379 }
5380 else
5381 {
5382 AssertDimension(n_data_sets, patches[0].data.n_rows())
5383 }
5384
5385 const char *ascii_or_binary =
5386 (deal_ii_with_zlib &&
5388 "binary" :
5389 "ascii";
5390
5391
5392 // first count the number of cells and cells for later use
5393 unsigned int n_nodes;
5394 unsigned int n_cells;
5395 std::tie(n_nodes, n_cells, std::ignore) =
5396 count_nodes_and_cells_and_points(patches, flags.write_higher_order_cells);
5397
5398 // -----------------
5399 // In the following, let us first set up a number of lambda functions that
5400 // will be used in building the different parts of the VTU file. We will
5401 // later call them in turn on different tasks.
5402 // first make up a list of used vertices along with their coordinates
5403 const auto stringize_vertex_information = [&patches,
5404 &flags,
5405 output_precision =
5406 out.precision(),
5407 ascii_or_binary]() {
5408 std::ostringstream o;
5409 o << " <Points>\n";
5410 o << " <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\""
5411 << ascii_or_binary << "\">\n";
5412 const std::vector<Point<spacedim>> node_positions =
5413 get_node_positions(patches);
5414
5415 // VTK/VTU always wants to see three coordinates, even if we are
5416 // in 1d or 2d. So pad node positions with zeros as appropriate.
5417 std::vector<float> node_coordinates_3d;
5418 node_coordinates_3d.reserve(node_positions.size() * 3);
5419 for (const auto &node_position : node_positions)
5420 {
5421 for (unsigned int d = 0; d < 3; ++d)
5422 if (d < spacedim)
5423 node_coordinates_3d.emplace_back(node_position[d]);
5424 else
5425 node_coordinates_3d.emplace_back(0.0f);
5426 }
5427 o << vtu_stringize_array(node_coordinates_3d,
5428 flags.compression_level,
5429 output_precision)
5430 << '\n';
5431 o << " </DataArray>\n";
5432 o << " </Points>\n\n";
5433
5434 return o.str();
5435 };
5436
5437
5438 //-------------------------------
5439 // Now for the cells. The first part of this is how vertices
5440 // build cells.
5441 const auto stringize_cell_to_vertex_information = [&patches,
5442 &flags,
5443 ascii_or_binary,
5444 output_precision =
5445 out.precision()]() {
5446 std::ostringstream o;
5447
5448 o << " <Cells>\n";
5449 o << " <DataArray type=\"Int32\" Name=\"connectivity\" format=\""
5450 << ascii_or_binary << "\">\n";
5451
5452 std::vector<int32_t> cells;
5453 Assert(dim <= 3, ExcNotImplemented());
5454
5455 unsigned int first_vertex_of_patch = 0;
5456
5457 for (const auto &patch : patches)
5458 {
5459 // First treat a slight oddball case: For triangles and tetrahedra,
5460 // the case with n_subdivisions==2 is treated as if the cell was
5461 // output as a single, quadratic, cell rather than as one would
5462 // expect as 4 sub-cells (for triangles; and the corresponding
5463 // number of sub-cells for tetrahedra). This is courtesy of some
5464 // special-casing in the function extract_vtk_patch_info().
5465 if ((dim >= 2) &&
5466 (patch.reference_cell == ReferenceCells::get_simplex<dim>()) &&
5467 (patch.n_subdivisions == 2))
5468 {
5469 const unsigned int n_points = patch.data.n_cols();
5470 Assert((dim == 2 && n_points == 6) ||
5471 (dim == 3 && n_points == 10),
5473
5474 if (deal_ii_with_zlib &&
5475 (flags.compression_level !=
5477 {
5478 for (unsigned int i = 0; i < n_points; ++i)
5479 cells.push_back(first_vertex_of_patch + i);
5480 }
5481 else
5482 {
5483 for (unsigned int i = 0; i < n_points; ++i)
5484 o << '\t' << first_vertex_of_patch + i;
5485 o << '\n';
5486 }
5487
5488 first_vertex_of_patch += n_points;
5489 }
5490 // Then treat all of the other non-hypercube cases since they can
5491 // currently not be subdivided (into sub-cells, or into higher-order
5492 // cells):
5493 else if (patch.reference_cell != ReferenceCells::get_hypercube<dim>())
5494 {
5496
5497 const unsigned int n_points = patch.data.n_cols();
5498
5499 if (deal_ii_with_zlib &&
5500 (flags.compression_level !=
5502 {
5503 for (unsigned int i = 0; i < n_points; ++i)
5504 cells.push_back(
5505 first_vertex_of_patch +
5507 }
5508 else
5509 {
5510 for (unsigned int i = 0; i < n_points; ++i)
5511 o << '\t'
5512 << (first_vertex_of_patch +
5514 o << '\n';
5515 }
5516
5517 first_vertex_of_patch += n_points;
5518 }
5519 else // a hypercube cell
5520 {
5521 const unsigned int n_subdivisions = patch.n_subdivisions;
5522 const unsigned int n_points_per_direction = n_subdivisions + 1;
5523
5524 std::vector<unsigned> local_vertex_order;
5525
5526 // Output the current state of the local_vertex_order array,
5527 // then clear it:
5528 const auto flush_current_cell = [&flags,
5529 &o,
5530 &cells,
5531 first_vertex_of_patch,
5532 &local_vertex_order]() {
5533 if (deal_ii_with_zlib &&
5534 (flags.compression_level !=
5536 {
5537 for (const auto &c : local_vertex_order)
5538 cells.push_back(first_vertex_of_patch + c);
5539 }
5540 else
5541 {
5542 for (const auto &c : local_vertex_order)
5543 o << '\t' << first_vertex_of_patch + c;
5544 o << '\n';
5545 }
5546
5547 local_vertex_order.clear();
5548 };
5549
5550 if (flags.write_higher_order_cells == false)
5551 {
5552 local_vertex_order.reserve(Utilities::fixed_power<dim>(2));
5553
5554 switch (dim)
5555 {
5556 case 0:
5557 {
5558 local_vertex_order.emplace_back(0);
5559 flush_current_cell();
5560 break;
5561 }
5562
5563 case 1:
5564 {
5565 for (unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
5566 {
5567 const unsigned int starting_offset = i1;
5568 local_vertex_order.emplace_back(starting_offset);
5569 local_vertex_order.emplace_back(starting_offset +
5570 1);
5571 flush_current_cell();
5572 }
5573 break;
5574 }
5575
5576 case 2:
5577 {
5578 for (unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
5579 for (unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
5580 {
5581 const unsigned int starting_offset =
5582 i2 * n_points_per_direction + i1;
5583 local_vertex_order.emplace_back(
5584 starting_offset);
5585 local_vertex_order.emplace_back(
5586 starting_offset + 1);
5587 local_vertex_order.emplace_back(
5588 starting_offset + n_points_per_direction + 1);
5589 local_vertex_order.emplace_back(
5590 starting_offset + n_points_per_direction);
5591 flush_current_cell();
5592 }
5593 break;
5594 }
5595
5596 case 3:
5597 {
5598 for (unsigned int i3 = 0; i3 < n_subdivisions; ++i3)
5599 for (unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
5600 for (unsigned int i1 = 0; i1 < n_subdivisions;
5601 ++i1)
5602 {
5603 const unsigned int starting_offset =
5604 i3 * n_points_per_direction *
5605 n_points_per_direction +
5606 i2 * n_points_per_direction + i1;
5607 local_vertex_order.emplace_back(
5608 starting_offset);
5609 local_vertex_order.emplace_back(
5610 starting_offset + 1);
5611 local_vertex_order.emplace_back(
5612 starting_offset + n_points_per_direction +
5613 1);
5614 local_vertex_order.emplace_back(
5615 starting_offset + n_points_per_direction);
5616 local_vertex_order.emplace_back(
5617 starting_offset + n_points_per_direction *
5618 n_points_per_direction);
5619 local_vertex_order.emplace_back(
5620 starting_offset +
5621 n_points_per_direction *
5622 n_points_per_direction +
5623 1);
5624 local_vertex_order.emplace_back(
5625 starting_offset +
5626 n_points_per_direction *
5627 n_points_per_direction +
5628 n_points_per_direction + 1);
5629 local_vertex_order.emplace_back(
5630 starting_offset +
5631 n_points_per_direction *
5632 n_points_per_direction +
5633 n_points_per_direction);
5634 flush_current_cell();
5635 }
5636 break;
5637 }
5638
5639 default:
5640 Assert(false, ExcNotImplemented());
5641 }
5642 }
5643 else // use higher-order output
5644 {
5645 local_vertex_order.resize(
5646 Utilities::fixed_power<dim>(n_points_per_direction));
5647
5648 switch (dim)
5649 {
5650 case 0:
5651 {
5652 Assert(false,
5653 ExcMessage(
5654 "Point-like cells should not be possible "
5655 "when writing higher-order cells."));
5656 break;
5657 }
5658 case 1:
5659 {
5660 for (unsigned int i1 = 0; i1 < n_subdivisions + 1;
5661 ++i1)
5662 {
5663 const unsigned int local_index = i1;
5664 const unsigned int connectivity_index =
5665 patch.reference_cell
5666 .template vtk_lexicographic_to_node_index<1>(
5667 {{i1}},
5668 {{n_subdivisions}},
5669 /* use VTU, not VTK: */ false);
5670 local_vertex_order[connectivity_index] =
5671 local_index;
5672 flush_current_cell();
5673 }
5674
5675 break;
5676 }
5677 case 2:
5678 {
5679 for (unsigned int i2 = 0; i2 < n_subdivisions + 1;
5680 ++i2)
5681 for (unsigned int i1 = 0; i1 < n_subdivisions + 1;
5682 ++i1)
5683 {
5684 const unsigned int local_index =
5685 i2 * n_points_per_direction + i1;
5686 const unsigned int connectivity_index =
5687 patch.reference_cell
5688 .template vtk_lexicographic_to_node_index<
5689 2>({{i1, i2}},
5690 {{n_subdivisions, n_subdivisions}},
5691 /* use VTU, not VTK: */ false);
5692 local_vertex_order[connectivity_index] =
5693 local_index;
5694 }
5695 flush_current_cell();
5696
5697 break;
5698 }
5699 case 3:
5700 {
5701 for (unsigned int i3 = 0; i3 < n_subdivisions + 1;
5702 ++i3)
5703 for (unsigned int i2 = 0; i2 < n_subdivisions + 1;
5704 ++i2)
5705 for (unsigned int i1 = 0; i1 < n_subdivisions + 1;
5706 ++i1)
5707 {
5708 const unsigned int local_index =
5709 i3 * n_points_per_direction *
5710 n_points_per_direction +
5711 i2 * n_points_per_direction + i1;
5712 const unsigned int connectivity_index =
5713 patch.reference_cell
5714 .template vtk_lexicographic_to_node_index<
5715 3>({{i1, i2, i3}},
5716 {{n_subdivisions,
5717 n_subdivisions,
5718 n_subdivisions}},
5719 /* use VTU, not VTK: */ false);
5720 local_vertex_order[connectivity_index] =
5721 local_index;
5722 }
5723
5724 flush_current_cell();
5725 break;
5726 }
5727 default:
5728 Assert(false, ExcNotImplemented());
5729 }
5730 }
5731
5732 // Finally update the number of the first vertex of this
5733 // patch
5734 first_vertex_of_patch +=
5735 Utilities::fixed_power<dim>(patch.n_subdivisions + 1);
5736 }
5737 }
5738
5739 // Flush the 'cells' object we created herein.
5740 if (deal_ii_with_zlib && (flags.compression_level !=
5742 {
5743 o << vtu_stringize_array(cells,
5744 flags.compression_level,
5745 output_precision)
5746 << '\n';
5747 }
5748 o << " </DataArray>\n";
5749
5750 return o.str();
5751 };
5752
5753
5754 //-------------------------------
5755 // The second part of cell information is the offsets in
5756 // the array built by the previous lambda function that indicate
5757 // individual cells.
5758 //
5759 // Note that this separates XML VTU format from the VTK format; the latter
5760 // puts the number of nodes per cell in front of the connectivity list for
5761 // each cell, whereas the VTU format uses one large list of vertex indices
5762 // and a separate array of offsets.
5763 //
5764 // The third piece to cell information is that we need to
5765 // output the types of the cells.
5766 //
5767 // The following function does both of these pieces.
5768 const auto stringize_cell_offset_and_type_information =
5769 [&patches,
5770 &flags,
5771 ascii_or_binary,
5772 n_cells,
5773 output_precision = out.precision()]() {
5774 std::ostringstream o;
5775
5776 o << " <DataArray type=\"Int32\" Name=\"offsets\" format=\""
5777 << ascii_or_binary << "\">\n";
5778
5779 std::vector<int32_t> offsets;
5780 offsets.reserve(n_cells);
5781
5782 // std::uint8_t might be an alias to unsigned char which is then not
5783 // printed as ascii integers
5784 std::vector<unsigned int> cell_types;
5785 cell_types.reserve(n_cells);
5786
5787 unsigned int first_vertex_of_patch = 0;
5788
5789 for (const auto &patch : patches)
5790 {
5791 const auto vtk_cell_id =
5792 extract_vtk_patch_info(patch, flags.write_higher_order_cells);
5793
5794 for (unsigned int i = 0; i < vtk_cell_id[1]; ++i)
5795 {
5796 cell_types.push_back(vtk_cell_id[0]);
5797 first_vertex_of_patch += vtk_cell_id[2];
5798 offsets.push_back(first_vertex_of_patch);
5799 }
5800 }
5801
5802 o << vtu_stringize_array(offsets,
5803 flags.compression_level,
5804 output_precision);
5805 o << '\n';
5806 o << " </DataArray>\n";
5807
5808 o << " <DataArray type=\"UInt8\" Name=\"types\" format=\""
5809 << ascii_or_binary << "\">\n";
5810
5811 if (deal_ii_with_zlib &&
5813 {
5814 std::vector<uint8_t> cell_types_uint8_t(cell_types.size());
5815 for (unsigned int i = 0; i < cell_types.size(); ++i)
5816 cell_types_uint8_t[i] = static_cast<std::uint8_t>(cell_types[i]);
5817
5818 o << vtu_stringize_array(cell_types_uint8_t,
5819 flags.compression_level,
5820 output_precision);
5821 }
5822 else
5823 {
5824 o << vtu_stringize_array(cell_types,
5825 flags.compression_level,
5826 output_precision);
5827 }
5828
5829 o << '\n';
5830 o << " </DataArray>\n";
5831 o << " </Cells>\n";
5832
5833 return o.str();
5834 };
5835
5836
5837 //-------------------------------------
5838 // data output.
5839
5840 const auto stringize_nonscalar_data_range =
5841 [&flags,
5842 &data_names,
5843 ascii_or_binary,
5844 n_data_sets,
5845 n_nodes,
5846 output_precision = out.precision()](const Table<2, float> &data_vectors,
5847 const auto & range) {
5848 std::ostringstream o;
5849
5850 const auto first_component = std::get<0>(range);
5851 const auto last_component = std::get<1>(range);
5852 const auto &name = std::get<2>(range);
5853 const bool is_tensor =
5854 (std::get<3>(range) ==
5856 const unsigned int n_components = (is_tensor ? 9 : 3);
5857 AssertThrow(last_component >= first_component,
5858 ExcLowerRange(last_component, first_component));
5859 AssertThrow(last_component < n_data_sets,
5860 ExcIndexRange(last_component, 0, n_data_sets));
5861 if (is_tensor)
5862 {
5863 AssertThrow((last_component + 1 - first_component <= 9),
5864 ExcMessage(
5865 "Can't declare a tensor with more than 9 components "
5866 "in VTK/VTU format."));
5867 }
5868 else
5869 {
5870 AssertThrow((last_component + 1 - first_component <= 3),
5871 ExcMessage(
5872 "Can't declare a vector with more than 3 components "
5873 "in VTK/VTU format."));
5874 }
5875
5876 // write the header. concatenate all the component names with double
5877 // underscores unless a vector name has been specified
5878 o << " <DataArray type=\"Float32\" Name=\"";
5879
5880 if (!name.empty())
5881 o << name;
5882 else
5883 {
5884 for (unsigned int i = first_component; i < last_component; ++i)
5885 o << data_names[i] << "__";
5886 o << data_names[last_component];
5887 }
5888
5889 o << "\" NumberOfComponents=\"" << n_components << "\" format=\""
5890 << ascii_or_binary << "\"";
5891 // If present, also list the physical units for this quantity. Look
5892 // this up for either the name of the whole vector/tensor, or if that
5893 // isn't listed, via its first component.
5894 if (!name.empty())
5895 {
5896 if (flags.physical_units.find(name) != flags.physical_units.end())
5897 o << " units=\"" << flags.physical_units.at(name) << "\"";
5898 }
5899 else
5900 {
5901 if (flags.physical_units.find(data_names[first_component]) !=
5902 flags.physical_units.end())
5903 o << " units=\""
5904 << flags.physical_units.at(data_names[first_component]) << "\"";
5905 }
5906 o << ">\n";
5907
5908 // now write data. pad all vectors to have three components
5909 std::vector<float> data;
5910 data.reserve(n_nodes * n_components);
5911
5912 for (unsigned int n = 0; n < n_nodes; ++n)
5913 {
5914 if (!is_tensor)
5915 {
5916 switch (last_component - first_component)
5917 {
5918 case 0:
5919 data.push_back(data_vectors(first_component, n));
5920 data.push_back(0);
5921 data.push_back(0);
5922 break;
5923
5924 case 1:
5925 data.push_back(data_vectors(first_component, n));
5926 data.push_back(data_vectors(first_component + 1, n));
5927 data.push_back(0);
5928 break;
5929
5930 case 2:
5931 data.push_back(data_vectors(first_component, n));
5932 data.push_back(data_vectors(first_component + 1, n));
5933 data.push_back(data_vectors(first_component + 2, n));
5934 break;
5935
5936 default:
5937 // Anything else is not yet implemented
5938 Assert(false, ExcInternalError());
5939 }
5940 }
5941 else
5942 {
5943 Tensor<2, 3> vtk_data;
5944 vtk_data = 0.;
5945
5946 const unsigned int size = last_component - first_component + 1;
5947 if (size == 1)
5948 // 1d, 1 element
5949 {
5950 vtk_data[0][0] = data_vectors(first_component, n);
5951 }
5952 else if (size == 4)
5953 // 2d, 4 elements
5954 {
5955 for (unsigned int c = 0; c < size; ++c)
5956 {
5957 const auto ind =
5959 vtk_data[ind[0]][ind[1]] =
5960 data_vectors(first_component + c, n);
5961 }
5962 }
5963 else if (size == 9)
5964 // 3d 9 elements
5965 {
5966 for (unsigned int c = 0; c < size; ++c)
5967 {
5968 const auto ind =
5970 vtk_data[ind[0]][ind[1]] =
5971 data_vectors(first_component + c, n);
5972 }
5973 }
5974 else
5975 {
5976 Assert(false, ExcInternalError());
5977 }
5978
5979 // now put the tensor into data
5980 // note we pad with zeros because VTK format always wants to
5981 // see a 3x3 tensor, regardless of dimension
5982 for (unsigned int i = 0; i < 3; ++i)
5983 for (unsigned int j = 0; j < 3; ++j)
5984 data.push_back(vtk_data[i][j]);
5985 }
5986 } // loop over nodes
5987
5988 o << vtu_stringize_array(data,
5989 flags.compression_level,
5990 output_precision);
5991 o << '\n';
5992 o << " </DataArray>\n";
5993
5994 return o.str();
5995 };
5996
5997 const auto stringize_scalar_data_set =
5998 [&flags,
5999 &data_names,
6000 ascii_or_binary,
6001 output_precision = out.precision()](const Table<2, float> &data_vectors,
6002 const unsigned int data_set) {
6003 std::ostringstream o;
6004
6005 o << " <DataArray type=\"Float32\" Name=\"" << data_names[data_set]
6006 << "\" format=\"" << ascii_or_binary << "\"";
6007 // If present, also list the physical units for this quantity.
6008 if (flags.physical_units.find(data_names[data_set]) !=
6009 flags.physical_units.end())
6010 o << " units=\"" << flags.physical_units.at(data_names[data_set])
6011 << "\"";
6012
6013 o << ">\n";
6014
6015 const std::vector<float> data(data_vectors[data_set].begin(),
6016 data_vectors[data_set].end());
6017 o << vtu_stringize_array(data,
6018 flags.compression_level,
6019 output_precision);
6020 o << '\n';
6021 o << " </DataArray>\n";
6022
6023 return o.str();
6024 };
6025
6026
6027 // For the format we write here, we need to write all node values relating
6028 // to one variable at a time. We could in principle do this by looping
6029 // over all patches and extracting the values corresponding to the one
6030 // variable we're dealing with right now, and then start the process over
6031 // for the next variable with another loop over all patches.
6032 //
6033 // An easier way is to create a global table that for each variable
6034 // lists all values. This copying of data vectors can be done in the
6035 // background while we're already working on vertices and cells,
6036 // so do this on a separate task and when wanting to write out the
6037 // data, we wait for that task to finish.
6039 create_global_data_table_task = Threads::new_task([&patches]() {
6040 return create_global_data_table<dim, spacedim, float>(patches);
6041 });
6042
6043 // -----------------------------
6044 // Now finally get around to actually doing anything. Let's start with
6045 // running the first three tasks generating the vertex and cell information:
6047 mesh_tasks += Threads::new_task(stringize_vertex_information);
6048 mesh_tasks += Threads::new_task(stringize_cell_to_vertex_information);
6049 mesh_tasks += Threads::new_task(stringize_cell_offset_and_type_information);
6050
6051 // For what follows, we have to have the reordered data available. So wait
6052 // for that task to conclude and get the resulting data table:
6053 const Table<2, float> data_vectors =
6054 std::move(*create_global_data_table_task.return_value());
6055
6056 // Then create the strings for the actual values of the solution vectors,
6057 // again on separate tasks:
6059 // When writing, first write out all vector and tensor data
6060 std::vector<bool> data_set_handled(n_data_sets, false);
6061 for (const auto &range : nonscalar_data_ranges)
6062 {
6063 // Mark these components as already handled:
6064 const auto first_component = std::get<0>(range);
6065 const auto last_component = std::get<1>(range);
6066 for (unsigned int i = first_component; i <= last_component; ++i)
6067 data_set_handled[i] = true;
6068
6069 data_tasks += Threads::new_task([&, range]() {
6070 return stringize_nonscalar_data_range(data_vectors, range);
6071 });
6072 }
6073
6074 // Now do the left over scalar data sets
6075 for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
6076 if (data_set_handled[data_set] == false)
6077 {
6078 data_tasks += Threads::new_task([&, data_set]() {
6079 return stringize_scalar_data_set(data_vectors, data_set);
6080 });
6081 }
6082
6083 // Alright, all tasks are now running. Wait for their conclusion and output
6084 // all of the data they have produced:
6085 out << "<Piece NumberOfPoints=\"" << n_nodes << "\" NumberOfCells=\""
6086 << n_cells << "\" >\n";
6087 for (const auto &s : mesh_tasks.return_values())
6088 out << s;
6089 out << " <PointData Scalars=\"scalars\">\n";
6090 for (const auto &s : data_tasks.return_values())
6091 out << s;
6092 out << " </PointData>\n";
6093 out << " </Piece>\n";
6094
6095 // make sure everything now gets to disk
6096 out.flush();
6097
6098 // assert the stream is still ok
6099 AssertThrow(out.fail() == false, ExcIO());
6100 }
6101
6102
6103
6104 void
6106 std::ostream & out,
6107 const std::vector<std::string> &piece_names,
6108 const std::vector<std::string> &data_names,
6109 const std::vector<
6110 std::tuple<unsigned int,
6111 unsigned int,
6112 std::string,
6114 & nonscalar_data_ranges,
6115 const VtkFlags &flags)
6116 {
6117 AssertThrow(out.fail() == false, ExcIO());
6118
6119 // If the user provided physical units, make sure that they don't contain
6120 // quote characters as this would make the VTU file invalid XML and
6121 // probably lead to all sorts of difficult error messages. Other than that,
6122 // trust the user that whatever they provide makes sense somehow.
6123 for (const auto &unit : flags.physical_units)
6124 {
6125 (void)unit;
6126 Assert(
6127 unit.second.find('\"') == std::string::npos,
6128 ExcMessage(
6129 "A physical unit you provided, <" + unit.second +
6130 ">, contained a quotation mark character. This is not allowed."));
6131 }
6132
6133 const unsigned int n_data_sets = data_names.size();
6134
6135 out << "<?xml version=\"1.0\"?>\n";
6136
6137 out << "<!--\n";
6138 out << "#This file was generated by the deal.II library"
6139 << " on " << Utilities::System::get_date() << " at "
6140 << Utilities::System::get_time() << "\n-->\n";
6141
6142 out
6143 << "<VTKFile type=\"PUnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
6144 out << " <PUnstructuredGrid GhostLevel=\"0\">\n";
6145 out << " <PPointData Scalars=\"scalars\">\n";
6146
6147 // We need to output in the same order as the write_vtu function does:
6148 std::vector<bool> data_set_written(n_data_sets, false);
6149 for (const auto &nonscalar_data_range : nonscalar_data_ranges)
6150 {
6151 const auto first_component = std::get<0>(nonscalar_data_range);
6152 const auto last_component = std::get<1>(nonscalar_data_range);
6153 const bool is_tensor =
6154 (std::get<3>(nonscalar_data_range) ==
6156 const unsigned int n_components = (is_tensor ? 9 : 3);
6157 AssertThrow(last_component >= first_component,
6158 ExcLowerRange(last_component, first_component));
6159 AssertThrow(last_component < n_data_sets,
6160 ExcIndexRange(last_component, 0, n_data_sets));
6161 if (is_tensor)
6162 {
6163 AssertThrow((last_component + 1 - first_component <= 9),
6164 ExcMessage(
6165 "Can't declare a tensor with more than 9 components "
6166 "in VTK"));
6167 }
6168 else
6169 {
6170 Assert((last_component + 1 - first_component <= 3),
6171 ExcMessage(
6172 "Can't declare a vector with more than 3 components "
6173 "in VTK"));
6174 }
6175
6176 // mark these components as already written:
6177 for (unsigned int i = std::get<0>(nonscalar_data_range);
6178 i <= std::get<1>(nonscalar_data_range);
6179 ++i)
6180 data_set_written[i] = true;
6181
6182 // write the header. concatenate all the component names with double
6183 // underscores unless a vector name has been specified
6184 out << " <PDataArray type=\"Float32\" Name=\"";
6185
6186 const std::string &name = std::get<2>(nonscalar_data_range);
6187 if (!name.empty())
6188 out << name;
6189 else
6190 {
6191 for (unsigned int i = std::get<0>(nonscalar_data_range);
6192 i < std::get<1>(nonscalar_data_range);
6193 ++i)
6194 out << data_names[i] << "__";
6195 out << data_names[std::get<1>(nonscalar_data_range)];
6196 }
6197
6198 out << "\" NumberOfComponents=\"" << n_components
6199 << "\" format=\"ascii\"";
6200 // If present, also list the physical units for this quantity. Look this
6201 // up for either the name of the whole vector/tensor, or if that isn't
6202 // listed, via its first component.
6203 if (!name.empty())
6204 {