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