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