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