1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2020 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 #ifndef dealii_data_out_dof_data_h
17 #define dealii_data_out_dof_data_h
18 
19 
20 
21 #include <deal.II/base/config.h>
22 
23 #include <deal.II/base/data_out_base.h>
24 #include <deal.II/base/mg_level_object.h>
25 #include <deal.II/base/smartpointer.h>
26 
27 #include <deal.II/dofs/dof_handler.h>
28 
29 #include <deal.II/fe/mapping.h>
30 
31 #include <deal.II/grid/tria.h>
32 
33 #include <deal.II/hp/fe_collection.h>
34 #include <deal.II/hp/fe_values.h>
35 #include <deal.II/hp/mapping_collection.h>
36 #include <deal.II/hp/q_collection.h>
37 
38 #include <deal.II/numerics/data_component_interpretation.h>
39 #include <deal.II/numerics/data_postprocessor.h>
40 
41 #include <memory>
42 
43 DEAL_II_NAMESPACE_OPEN
44 
45 namespace Exceptions
46 {
47   /**
48    * A namespace for exceptions that are used throughout the DataOut*
49    * collection of classes.
50    */
51   namespace DataOutImplementation
52   {
53     /**
54      * Exception
55      */
56     DeclException1(ExcInvalidNumberOfSubdivisions,
57                    int,
58                    << "The number of subdivisions per patch, " << arg1
59                    << ", is not valid. It needs to be greater or equal to "
60                       "one, or zero if you want it to be determined "
61                       "automatically.");
62 
63     /**
64      * Exception
65      */
66     DeclExceptionMsg(ExcNoTriangulationSelected,
67                      "For the operation you are attempting, you first need to "
68                      "tell the DataOut or related object which DoFHandler or "
69                      "triangulation you would like to work on.");
70 
71     /**
72      * Exception
73      */
74     DeclExceptionMsg(ExcNoDoFHandlerSelected,
75                      "For the operation you are attempting, you first need to "
76                      "tell the DataOut or related object which DoFHandler "
77                      "you would like to work on.");
78 
79     /**
80      * Exception
81      */
82     DeclException3(ExcInvalidVectorSize,
83                    int,
84                    int,
85                    int,
86                    << "The vector has size " << arg1
87                    << " but the DoFHandler object says that there are " << arg2
88                    << " degrees of freedom and there are " << arg3
89                    << " active cells. The size of your vector needs to be"
90                    << " either equal to the number of degrees of freedom (when"
91                    << " the data is of type type_dof_data), or equal to the"
92                    << " number of active cells (when the data is of type "
93                    << " type_cell_data).");
94     /**
95      * Exception
96      */
97     DeclException2(
98       ExcInvalidCharacter,
99       std::string,
100       size_t,
101       << "Please use only the characters [a-zA-Z0-9_<>()] for" << std::endl
102       << "description strings since some graphics formats will only accept these."
103       << std::endl
104       << "The string you gave was <" << arg1
105       << ">, within which the invalid character is <" << arg1[arg2] << ">."
106       << std::endl);
107     /**
108      * Exception
109      */
110     DeclExceptionMsg(
111       ExcOldDataStillPresent,
112       "When attaching a triangulation or DoFHandler object, it is "
113       "not allowed if old data vectors are still referenced. If "
114       "you want to reuse an object of the current type, you first "
115       "need to call the 'clear_data_vector()' function.");
116     /**
117      * Exception
118      */
119     DeclException2(ExcInvalidNumberOfNames,
120                    int,
121                    int,
122                    << "You have to give one name per component in your "
123                    << "data vector. The number you gave was " << arg1
124                    << ", but the number of components is " << arg2 << ".");
125     /**
126      * Exception
127      */
128     DeclExceptionMsg(ExcIncompatibleDatasetNames,
129                      "While merging sets of patches, the two sets to be merged "
130                      "need to refer to data that agrees on the names of the "
131                      "various variables represented. In other words, you "
132                      "cannot merge sets of patches that originate from "
133                      "entirely unrelated simulations.");
134     /**
135      * Exception
136      */
137     DeclExceptionMsg(ExcIncompatiblePatchLists,
138                      "While merging sets of patches, the two sets to be merged "
139                      "need to refer to data that agrees on the number of "
140                      "subdivisions and other properties. In other words, you "
141                      "cannot merge sets of patches that originate from "
142                      "entirely unrelated simulations.");
143 
144     DeclException2(ExcInvalidVectorDeclaration,
145                    int,
146                    std::string,
147                    << "When declaring that a number of components in a data "
148                    << "set to be output logically form a vector instead of "
149                    << "simply a set of scalar fields, you need to specify "
150                    << "this for all relevant components. Furthermore, "
151                    << "vectors must always consist of exactly <dim> "
152                    << "components. However, the vector component at "
153                    << "position " << arg1 << " with name <" << arg2
154                    << "> does not satisfy these conditions.");
155 
156     DeclException2(ExcInvalidTensorDeclaration,
157                    int,
158                    std::string,
159                    << "When declaring that a number of components in a data "
160                    << "set to be output logically form a tensor instead of "
161                    << "simply a set of scalar fields, you need to specify "
162                    << "this for all relevant components. Furthermore, "
163                    << "tensors must always consist of exactly <dim*dim> "
164                    << "components. However, the tensor component at "
165                    << "position " << arg1 << " with name <" << arg2
166                    << "> does not satisfy these conditions.");
167 
168   } // namespace DataOutImplementation
169 } // namespace Exceptions
170 
171 
172 namespace internal
173 {
174   namespace DataOutImplementation
175   {
176     /**
177      * The DataEntry classes abstract away the concrete data type of vectors
178      * users can attach to DataOut (and similar) objects and allow the
179      * underlying DataOut functions to query for individual elements of solution
180      * vectors without having to know the concrete vector type. This avoids that
181      * DataOut has to know what vectors are being used, but it has the downside
182      * that DataOut also doesn't know the underlying scalar type of these
183      * vectors.
184      *
185      * If the underlying scalar types all represent real numbers (in the
186      * mathematical sense -- i.e., the scalar type would be @p float,
187      * @p double, etc) then that is not a problem -- DataOut simply
188      * receives the values of individual vector components as @p double
189      * objects. On the other hand, if the vector type uses a std::complex
190      * scalar type, then DataEntry returning a @p double for a vector
191      * entry is not sufficient -- we need to provide DataOut with a way
192      * to query both the real and the imaginary part, so that they can
193      * be written into output files separately.
194      *
195      * This enum allows DataOut to tell a DataEntry function which component
196      * of a vector entry it wants to query, i.e., whether it wants the real
197      * or the imaginary part of a vector entry.
198      */
199     enum class ComponentExtractor
200     {
201       real_part,
202       imaginary_part
203     };
204 
205 
206     /**
207      * For each vector that has been added through the add_data_vector()
208      * functions, we need to keep track of a pointer to it, and allow data
209      * extraction from it when we generate patches. Unfortunately, we need to
210      * do this for a number of different vector types. Fortunately, they all
211      * have the same interface. So the way we go is to have a base class that
212      * provides the functions to access the vector's information, and to have
213      * a derived template class that can be instantiated for each vector type.
214      * Since the vectors all have the same interface, this is no big problem,
215      * as they can all use the same general templatized code.
216      *
217      * @note This class is an example of the
218      * <a href="https://www.artima.com/cppsource/type_erasure.html">type
219      * erasure</a> design pattern.
220      */
221     template <typename DoFHandlerType>
222     class DataEntryBase
223     {
224     public:
225       /**
226        * Constructor. Give a list of names for the individual components of
227        * the vector and their interpretation as scalar or vector data. This
228        * constructor assumes that no postprocessor is going to be used.
229        */
230       DataEntryBase(const DoFHandlerType *          dofs,
231                     const std::vector<std::string> &names,
232                     const std::vector<
233                       DataComponentInterpretation::DataComponentInterpretation>
234                       &data_component_interpretation);
235 
236       /**
237        * Constructor when a data postprocessor is going to be used. In that
238        * case, the names and vector declarations are going to be acquired from
239        * the postprocessor.
240        */
241       DataEntryBase(const DoFHandlerType *dofs,
242                     const DataPostprocessor<DoFHandlerType::space_dimension>
243                       *data_postprocessor);
244 
245       /**
246        * Destructor made virtual.
247        */
248       virtual ~DataEntryBase() = default;
249 
250       /**
251        * Assuming that the stored vector is a cell vector, extract the given
252        * element from it.
253        */
254       virtual double
255       get_cell_data_value(const unsigned int       cell_number,
256                           const ComponentExtractor extract_component) const = 0;
257 
258       /**
259        * Given a FEValuesBase object, extract the values on the present cell
260        * from the vector we actually store.
261        */
262       virtual void
263       get_function_values(
264         const FEValuesBase<DoFHandlerType::dimension,
265                            DoFHandlerType::space_dimension> &fe_patch_values,
266         const ComponentExtractor                             extract_component,
267         std::vector<double> &patch_values) const = 0;
268 
269       /**
270        * Given a FEValuesBase object, extract the values on the present cell
271        * from the vector we actually store. This function does the same as the
272        * one above but for vector-valued finite elements.
273        */
274       virtual void
275       get_function_values(
276         const FEValuesBase<DoFHandlerType::dimension,
277                            DoFHandlerType::space_dimension> &fe_patch_values,
278         const ComponentExtractor                             extract_component,
279         std::vector<dealii::Vector<double>> &patch_values_system) const = 0;
280 
281       /**
282        * Given a FEValuesBase object, extract the gradients on the present
283        * cell from the vector we actually store.
284        */
285       virtual void
286       get_function_gradients(
287         const FEValuesBase<DoFHandlerType::dimension,
288                            DoFHandlerType::space_dimension> &fe_patch_values,
289         const ComponentExtractor                             extract_component,
290         std::vector<Tensor<1, DoFHandlerType::space_dimension>>
291           &patch_gradients) const = 0;
292 
293       /**
294        * Given a FEValuesBase object, extract the gradients on the present
295        * cell from the vector we actually store. This function does the same
296        * as the one above but for vector-valued finite elements.
297        */
298       virtual void
299       get_function_gradients(
300         const FEValuesBase<DoFHandlerType::dimension,
301                            DoFHandlerType::space_dimension> &fe_patch_values,
302         const ComponentExtractor                             extract_component,
303         std::vector<std::vector<Tensor<1, DoFHandlerType::space_dimension>>>
304           &patch_gradients_system) const = 0;
305 
306       /**
307        * Given a FEValuesBase object, extract the second derivatives on the
308        * present cell from the vector we actually store.
309        */
310       virtual void
311       get_function_hessians(
312         const FEValuesBase<DoFHandlerType::dimension,
313                            DoFHandlerType::space_dimension> &fe_patch_values,
314         const ComponentExtractor                             extract_component,
315         std::vector<Tensor<2, DoFHandlerType::space_dimension>> &patch_hessians)
316         const = 0;
317 
318       /**
319        * Given a FEValuesBase object, extract the second derivatives on the
320        * present cell from the vector we actually store. This function does
321        * the same as the one above but for vector-valued finite elements.
322        */
323       virtual void
324       get_function_hessians(
325         const FEValuesBase<DoFHandlerType::dimension,
326                            DoFHandlerType::space_dimension> &fe_patch_values,
327         const ComponentExtractor                             extract_component,
328         std::vector<std::vector<Tensor<2, DoFHandlerType::space_dimension>>>
329           &patch_hessians_system) const = 0;
330 
331       /**
332        * Return whether the data represented by (a derived class of) this object
333        * represents a complex-valued (as opposed to real-valued) information.
334        */
335       virtual bool
336       is_complex_valued() const = 0;
337 
338       /**
339        * Clear all references to the vectors.
340        */
341       virtual void
342       clear() = 0;
343 
344       /**
345        * Determine an estimate for the memory consumption (in bytes) of this
346        * object.
347        */
348       virtual std::size_t
349       memory_consumption() const = 0;
350 
351       /**
352        * Pointer to the DoFHandler object that the vector is based on.
353        */
354       SmartPointer<const DoFHandlerType> dof_handler;
355 
356       /**
357        * Names of the components of this data vector.
358        */
359       const std::vector<std::string> names;
360 
361       /**
362        * A vector that for each of the n_output_variables variables of the
363        * current data set indicates whether they are scalar fields, parts of a
364        * vector-field, or any of the other supported kinds of data.
365        */
366       const std::vector<
367         DataComponentInterpretation::DataComponentInterpretation>
368         data_component_interpretation;
369 
370       /**
371        * Pointer to a DataPostprocessing object which shall be applied to this
372        * data vector.
373        */
374       SmartPointer<
375         const dealii::DataPostprocessor<DoFHandlerType::space_dimension>>
376         postprocessor;
377 
378       /**
379        * Number of output variables this dataset provides (either number of
380        * components in vector valued function / data vector or number of
381        * computed quantities, if DataPostprocessor is applied). This variable
382        * is determined via and thus equivalent to <tt>names.size()</tt>.
383        */
384       unsigned int n_output_variables;
385     };
386 
387 
388     /**
389      * A data structure that holds all data needed in one thread when building
390      * patches in parallel. These data structures are created globally rather
391      * than on each cell to avoid allocation of memory in the threads. This is
392      * a base class for the AdditionalData kind of data structure discussed in
393      * the documentation of the WorkStream class.
394      *
395      * The <code>cell_to_patch_index_map</code> is an array that stores for
396      * index <tt>[i][j]</tt> the number of the patch that associated with the
397      * cell with index @p j on level @p i. This information is set up prior to
398      * generation of the patches, and is needed to generate neighborship
399      * information.
400      *
401      * This structure is used by several of the DataOut* classes, which
402      * derived their own ParallelData classes from it for additional fields.
403      */
404     template <int dim, int spacedim>
405     struct ParallelDataBase
406     {
407       ParallelDataBase(
408         const unsigned int               n_datasets,
409         const unsigned int               n_subdivisions,
410         const std::vector<unsigned int> &n_postprocessor_outputs,
411         const Mapping<dim, spacedim> &   mapping,
412         const std::vector<
413           std::shared_ptr<dealii::hp::FECollection<dim, spacedim>>>
414           &               finite_elements,
415         const UpdateFlags update_flags,
416         const bool        use_face_values);
417 
418       ParallelDataBase(
419         const unsigned int               n_datasets,
420         const unsigned int               n_subdivisions,
421         const std::vector<unsigned int> &n_postprocessor_outputs,
422         const dealii::hp::MappingCollection<dim, spacedim> &mapping,
423         const std::vector<
424           std::shared_ptr<dealii::hp::FECollection<dim, spacedim>>>
425           &               finite_elements,
426         const UpdateFlags update_flags,
427         const bool        use_face_values);
428 
429       ParallelDataBase(const ParallelDataBase &data);
430 
431       template <typename DoFHandlerType>
432       void
433       reinit_all_fe_values(
434         std::vector<std::shared_ptr<DataEntryBase<DoFHandlerType>>> &dof_data,
435         const typename dealii::Triangulation<dim, spacedim>::cell_iterator
436           &                cell,
437         const unsigned int face = numbers::invalid_unsigned_int);
438 
439       const FEValuesBase<dim, spacedim> &
440       get_present_fe_values(const unsigned int dataset) const;
441 
442       void
443       resize_system_vectors(const unsigned int n_components);
444 
445       const unsigned int n_datasets;
446       const unsigned int n_subdivisions;
447 
448       DataPostprocessorInputs::Scalar<spacedim>        patch_values_scalar;
449       DataPostprocessorInputs::Vector<spacedim>        patch_values_system;
450       std::vector<std::vector<dealii::Vector<double>>> postprocessed_values;
451 
452       const dealii::hp::MappingCollection<dim, spacedim> mapping_collection;
453       const std::vector<
454         std::shared_ptr<dealii::hp::FECollection<dim, spacedim>>>
455                         finite_elements;
456       const UpdateFlags update_flags;
457 
458       std::vector<std::shared_ptr<dealii::hp::FEValues<dim, spacedim>>>
459         x_fe_values;
460       std::vector<std::shared_ptr<dealii::hp::FEFaceValues<dim, spacedim>>>
461         x_fe_face_values;
462     };
463   } // namespace DataOutImplementation
464 } // namespace internal
465 
466 
467 // TODO: Most of the documentation of DataOut_DoFData applies to DataOut.
468 
469 /**
470  * This is an abstract class which provides the functionality to generate
471  * patches for output by base classes from data vectors on a grid. It allows
472  * to attach one or more pointers to a DoFHandler and attached node and cell
473  * data denoting functions on the grid which shall later be written in any of
474  * the implemented data formats.
475  *
476  *
477  * <h3>User visible interface</h3>
478  *
479  * The user visible interface of this class allows the user to specify data in
480  * two different ways. One is to make a DoFHandler object known to this class
481  * and to add data vectors that all correspond to this DoFHandler or the grid
482  * cells which will later be written to a file in some format. The second
483  * approach is to pass a DoFHandler object along with the vector. This allows
484  * setting data from different DoFHandlers in a neat way (of course, they both
485  * need to be based on the same triangulation). Instead of pondering about the
486  * different functions, an example for the first kind is probably the best
487  * explanation:
488  * @code
489  *   ...
490  *   ...   // compute solution, which contains nodal values
491  *   ...
492  *   ...   // compute error_estimator, which contains one value per cell
493  *
494  *   std::vector<std::string> solution_names;
495  *   solution_names.emplace_back ("x-displacement");
496  *   solution_names.emplace_back ("y-displacement");
497  *
498  *   DataOut<dim> data_out;
499  *   data_out.attach_dof_handler (dof_handler);
500  *   data_out.add_data_vector (solution, solution_names);
501  *   data_out.add_data_vector (error_estimator, "estimated_error");
502  *
503  *   data_out.build_patches ();
504  *
505  *   ofstream output_file ("output");
506  *   data_out.write_xxx (output_file);
507  *
508  *   data_out.clear();
509  * @endcode
510  *
511  * attach_dof_handler() tells this class that all future operations are to
512  * take place with the DoFHandler object and the triangulation it lives on. We
513  * then add the solution vector and the error estimator; note that they have
514  * different dimensions, because the solution is a nodal vector, here
515  * consisting of two components ("x-displacement" and "y-displacement") while
516  * the error estimator probably is a vector holding cell data. When attaching
517  * a data vector, you have to give a name to each component of the vector,
518  * which is done through an object of type <tt>vector<string></tt> as second
519  * argument; if only one component is in the vector, for example if we are
520  * adding cell data as in the second case, or if the finite element used by
521  * the DoFHandler has only one component, then you can use the second
522  * add_data_vector() function which takes a @p string instead of the
523  * <tt>vector<string></tt>.
524  *
525  * The add_data_vector() functions have additional arguments (with default
526  * values) that can be used to specify certain transformations. In particular,
527  * it allows to attach DataPostprocessor arguments to compute derived
528  * information from a data vector at each point at which the field will be
529  * evaluated so that it can be written to a file (for example, the Mach number
530  * in hypersonic flow can be computed from density and velocities; step-29
531  * also shows an example); another piece of information specified through
532  * arguments with default values is how certain output components should be
533  * interpreted, i.e. whether each component of the data is logically an
534  * independent scalar field, or whether some of them together form logically a
535  * vector-field (see the
536  * DataComponentInterpretation::DataComponentInterpretation enum, and the
537  * @ref step_22 "step-22"
538  * tutorial program).
539  *
540  * This class does not copy the vector given to it through the
541  * add_data_vector() functions, for memory consumption reasons. It only stores
542  * a reference to it, so it is in your responsibility to make sure that the
543  * data vectors exist long enough.
544  *
545  * After adding all data vectors, you need to call a function which generates
546  * the patches (i.e., some intermediate data representation) for output from
547  * the stored data. Derived classes name this function build_patches().
548  * Finally, you write() the data in one format or other, to a file.
549  *
550  * In the example above, an object of type DataOut was used, i.e. an object of
551  * a derived class. This is necessary since the current class does not provide
552  * means to actually generate the patches, only aids to store and access data.
553  * Any real functionality is implemented in derived classes such as DataOut.
554  *
555  * Note that the base class of this class, DataOutInterface offers several
556  * functions to ease programming with run-time determinable output formats
557  * (i.e. you need not use a fixed format by calling
558  * DataOutInterface::write_xxx in the above example, but you can select it by
559  * a run-time parameter without having to write the <tt>if () ... else
560  * ...</tt> clauses yourself), and also functions and classes offering ways to
561  * control the appearance of the output by setting flags for each output
562  * format.
563  *
564  *
565  * <h3>Information for derived classes</h3>
566  *
567  * What this class lacks is a way to produce the patches for output itself,
568  * from the stored data and degree of freedom information. Since this task is
569  * often application dependent it is left to derived classes. For example, in
570  * many applications, it might be wanted to limit the depth of output to a
571  * certain number of refinement levels and write data from finer cells only in
572  * a way interpolated to coarser cells, to reduce the amount of output. Also,
573  * it might be wanted to use different numbers of subdivisions on different
574  * cells when forming a patch, for example to accomplish for different
575  * polynomial degrees of the trial space on different cells. Also, the output
576  * need not necessarily consist of a patch for each cell, but might be made up
577  * of patches for faces, of other things. Take a look at derived classes to
578  * what is possible in this respect.
579  *
580  * For this reason, it is left to a derived class to provide a function, named
581  * usually build_patches() or the like, which fills the #patches array of this
582  * class.
583  *
584  * Regarding the templates of this class, it needs three values: first the
585  * space dimension in which the triangulation and the DoF handler operate,
586  * second the dimension of the objects which the patches represent.  Although
587  * in most cases they are equal, there are also classes for which this does
588  * not hold, for example if one outputs the result of a computation exploiting
589  * rotational symmetry in the original domain (in which the space dimension of
590  * the output would be one higher than that of the DoF handler, see the
591  * DataOut_Rotation() class), or one might conceive that one could write a
592  * class that only outputs the solution on a cut through the domain, in which
593  * case the space dimension of the output is less than that of the DoF
594  * handler. The last template argument denotes the dimension of the space into
595  * which the patches are embedded; usually, this dimension is the same as the
596  * dimensio of the patches themselves (which is also the default value of the
597  * template parameter), but there might be cases where this is not so. For
598  * example, in the DataOut_Faces() class, patches are generated from faces of
599  * the triangulation. Thus, the dimension of the patch is one less than the
600  * dimension of the embedding space, which is, in this case, equal to the
601  * dimension of the triangulation and DoF handler. However, for the cut
602  * through the domain mentioned above, if the cut is a straight one, then the
603  * cut can be embedded into a space of one dimension lower than the dimension
604  * of the triangulation, so that the last template parameter has the same
605  * value as the second one.
606  *
607  * @ingroup output
608  */
609 template <typename DoFHandlerType,
610           int patch_dim,
611           int patch_space_dim = patch_dim>
612 class DataOut_DoFData : public DataOutInterface<patch_dim, patch_space_dim>
613 {
614 public:
615   /**
616    * Typedef to the iterator type of the dof handler class under
617    * consideration.
618    */
619   using cell_iterator =
620     typename Triangulation<DoFHandlerType::dimension,
621                            DoFHandlerType::space_dimension>::cell_iterator;
622 
623 public:
624   /**
625    * Type describing what the vector given to add_data_vector() is: a vector
626    * that has one entry per degree of freedom in a DoFHandler object (such as
627    * solution vectors), or one entry per cell in the triangulation underlying
628    * the DoFHandler object (such as error per cell data). The value
629    * #type_automatic tells add_data_vector() to find out itself (see the
630    * documentation of add_data_vector() for the method used).
631    */
632   enum DataVectorType
633   {
634     /**
635      * Data vector entries are associated to degrees of freedom
636      */
637     type_dof_data,
638 
639     /**
640      * Data vector entries are one per grid cell
641      */
642     type_cell_data,
643 
644     /**
645      * Find out automatically
646      */
647     type_automatic
648   };
649 
650   /**
651    * Constructor
652    */
653   DataOut_DoFData();
654 
655   /**
656    * Destructor.
657    */
658   virtual ~DataOut_DoFData() override;
659 
660   /**
661    * Designate a dof handler to be used to extract geometry data and the
662    * mapping between nodes and node values. This call is not necessary if all
663    * added data vectors are supplemented with a DoFHandler argument.
664    *
665    * This call is optional: If you add data vectors with specified DoFHandler
666    * object, then that contains all information needed to generate the output.
667    */
668   void
669   attach_dof_handler(const DoFHandlerType &);
670 
671   /**
672    * Designate a triangulation to be used to extract geometry data and the
673    * mapping between nodes and node values.
674    *
675    * This call is optional: If you add data vectors with specified DoFHandler
676    * object, then that contains all information needed to generate the output.
677    * This call is useful when you only output cell vectors and no DoFHandler
678    * at all, in which case it provides the geometry.
679    */
680   void
681   attach_triangulation(const Triangulation<DoFHandlerType::dimension,
682                                            DoFHandlerType::space_dimension> &);
683 
684   /**
685    * Add a data vector together with its name.
686    *
687    * A pointer to the vector is stored, so you have to make sure the vector
688    * exists at that address at least as long as you call the <tt>write_*</tt>
689    * functions.
690    *
691    * It is assumed that the vector has the same number of components as there
692    * are degrees of freedom in the dof handler, in which case it is assumed to
693    * be a vector storing nodal data; or the size may be the number of active
694    * cells on the present grid, in which case it is assumed to be a cell data
695    * vector. As the number of degrees of freedom and of cells is usually not
696    * equal, the function can determine itself which type of vector it is
697    * given. However, there are corner cases where this automatic determination
698    * does not work.  One example is if you compute with piecewise constant
699    * elements and have a scalar solution, then there are as many cells as
700    * there are degrees of freedom (though they may be numbered differently).
701    * Another possibility is if you have a 1d mesh embedded in 2d space and the
702    * mesh consists of a closed curve of cells; in this case, there are as many
703    * nodes as there are cells, and when using a Q1 element you will have as
704    * many degrees of freedom as there are cells.  In these cases, you can
705    * change the last argument of the function from its default value
706    * #type_automatic to either #type_dof_data or #type_cell_data, depending on
707    * what the vector represents. Apart from such corner cases, you can leave
708    * the argument at its default value and let the function determine the type
709    * of the vector itself.
710    *
711    * If it is a vector holding DoF data, the names given shall be one for each
712    * component of the underlying finite element.  If it is a finite element
713    * composed of only one subelement, then there is another function following
714    * which takes a single name instead of a vector of names.
715    *
716    * The data_component_interpretation argument contains information about how
717    * the individual components of output files that consist of more than one
718    * data set are to be interpreted.
719    *
720    * For example, if one has a finite element for the Stokes equations in 2d,
721    * representing components (u,v,p), one would like to indicate that the
722    * first two, u and v, represent a logical vector so that later on when we
723    * generate graphical output we can hand them off to a visualization program
724    * that will automatically know to render them as a vector field, rather
725    * than as two separate and independent scalar fields.
726    *
727    * The default value of this argument (i.e. an empty vector) corresponds is
728    * equivalent to a vector of values
729    * DataComponentInterpretation::component_is_scalar, indicating that all
730    * output components are independent scalar fields. However, if the given
731    * data vector represents logical vectors, you may pass a vector that
732    * contains values DataComponentInterpretation::component_is_part_of_vector.
733    * In the example above, one would pass in a vector with components
734    * (DataComponentInterpretation::component_is_part_of_vector,
735    * DataComponentInterpretation::component_is_part_of_vector,
736    * DataComponentInterpretation::component_is_scalar) for (u,v,p).
737    *
738    * The names of a data vector shall only contain characters which are
739    * letters, underscore and a few other ones. Refer to the
740    * ExcInvalidCharacter exception declared in this class to see which
741    * characters are valid and which are not.
742    *
743    * @note The actual type for the vector argument may be any vector type from
744    * which FEValues can extract values on a cell using the
745    * FEValuesBase::get_function_values() function.
746    *
747    * @note When working in parallel, the vector to be written needs to be ghosted
748    * with read access to all degrees of freedom on the locally owned cells, see
749    * the step-40 or step-37 tutorial programs for details, i.e., it might be
750    * necessary to call data.update_ghost_values().
751    */
752   template <class VectorType>
753   void
754   add_data_vector(
755     const VectorType &              data,
756     const std::vector<std::string> &names,
757     const DataVectorType            type = type_automatic,
758     const std::vector<DataComponentInterpretation::DataComponentInterpretation>
759       &data_component_interpretation = std::vector<
760         DataComponentInterpretation::DataComponentInterpretation>());
761 
762   /**
763    * This function is an abbreviation to the above one (see there for a
764    * discussion of the various arguments), intended for use with finite
765    * elements that are not composed of subelements. In this case, only one
766    * name per data vector needs to be given, which is what this function
767    * takes. It simply relays its arguments after a conversion of the @p name
768    * to a vector of strings, to the other add_data_vector() function above.
769    *
770    * If @p data is a vector with multiple components this function will
771    * generate distinct names for all components by appending an underscore and
772    * the number of each component to @p name
773    *
774    * The actual type for the template argument may be any vector type from
775    * which FEValues can extract values on a cell using the
776    * FEValuesBase::get_function_values() function.
777    */
778   template <class VectorType>
779   void
780   add_data_vector(
781     const VectorType &   data,
782     const std::string &  name,
783     const DataVectorType type = type_automatic,
784     const std::vector<DataComponentInterpretation::DataComponentInterpretation>
785       &data_component_interpretation = std::vector<
786         DataComponentInterpretation::DataComponentInterpretation>());
787 
788   /**
789    * This function is an extension of the above one (see there for a
790    * discussion of the arguments except the first one) and allows to set a
791    * vector with its own DoFHandler object. This DoFHandler needs to be
792    * compatible with the other DoFHandler objects assigned with calls to @p
793    * add_data_vector or @p attach_dof_handler, in the sense that all of the
794    * DoFHandler objects need to be based on the same triangulation. This
795    * function allows you to export data from multiple DoFHandler objects that
796    * describe different solution components. An example of using this function
797    * is given in step-61.
798    *
799    * Since this function takes a DoFHandler object and hence naturally
800    * represents dof data, the data vector type argument present in the other
801    * methods above is not necessary.
802    */
803   template <class VectorType>
804   void
805   add_data_vector(
806     const DoFHandlerType &          dof_handler,
807     const VectorType &              data,
808     const std::vector<std::string> &names,
809     const std::vector<DataComponentInterpretation::DataComponentInterpretation>
810       &data_component_interpretation = std::vector<
811         DataComponentInterpretation::DataComponentInterpretation>());
812 
813 
814   /**
815    * This function is an abbreviation of the function above with only a scalar
816    * @p dof_handler given and a single data name.
817    */
818   template <class VectorType>
819   void
820   add_data_vector(
821     const DoFHandlerType &dof_handler,
822     const VectorType &    data,
823     const std::string &   name,
824     const std::vector<DataComponentInterpretation::DataComponentInterpretation>
825       &data_component_interpretation = std::vector<
826         DataComponentInterpretation::DataComponentInterpretation>());
827 
828   /**
829    * This function is an alternative to the above ones, allowing the output of
830    * derived quantities instead of the given data. This conversion has to be
831    * done in a class derived from DataPostprocessor. This function is used in
832    * step-29. Other uses are shown in step-32 and step-33.
833    *
834    * The names for these derived quantities are provided by the @p
835    * data_postprocessor argument. Likewise, the data_component_interpretation
836    * argument of the other add_data_vector() functions is provided by the
837    * data_postprocessor argument. As only data of type @p type_dof_data can be
838    * transformed, this type is also known implicitly and does not have to be
839    * given.
840    *
841    * @note The actual type for the vector argument may be any vector type from
842    * which FEValues can extract values on a cell using the
843    * FEValuesBase::get_function_values() function.
844    *
845    * @note The DataPostprocessor object (i.e., in reality the object of your
846    * derived class) has to live until the DataOut object is destroyed as the
847    * latter keeps a pointer to the former and will complain if the object
848    * pointed to is destroyed while the latter still has a pointer to it. If
849    * both the data postprocessor and DataOut objects are local variables of a
850    * function (as they are, for example, in step-29), then you can avoid this
851    * error by declaring the data postprocessor variable before the DataOut
852    * variable as objects are destroyed in reverse order of declaration.
853    */
854   template <class VectorType>
855   void
856   add_data_vector(const VectorType &data,
857                   const DataPostprocessor<DoFHandlerType::space_dimension>
858                     &data_postprocessor);
859 
860   /**
861    * Same function as above, but with a DoFHandler object that does not need
862    * to coincide with the DoFHandler initially set. Note that the
863    * postprocessor can only read data from the given DoFHandler and solution
864    * vector, not other solution vectors or DoFHandlers.
865    */
866   template <class VectorType>
867   void
868   add_data_vector(const DoFHandlerType &dof_handler,
869                   const VectorType &    data,
870                   const DataPostprocessor<DoFHandlerType::space_dimension>
871                     &data_postprocessor);
872 
873   /**
874    * Add a multilevel data vector.
875    *
876    * This function adds the vector-valued multilevel vector @p data in the
877    * form of a vector on each level that belongs to the DoFHandler @p
878    * dof_handler to the graphical output. This function is typically used in
879    * conjunction with a call to set_cell_selection() that selects cells on a
880    * specific level and not the active cells (the default).
881    *
882    * A vector @p data can be obtained in several ways, for example by using
883    * Multigrid::solution or Multigrid::defect during or after a multigrid
884    * cycle or by interpolating a solution via
885    * MGTransferMatrixFree::interpolate_to_mg().
886    *
887    * The handling of @p names and @p data_component_interpretation is identical
888    * to the add_data_vector() function.
889    */
890   template <class VectorType>
891   void
892   add_mg_data_vector(
893     const DoFHandlerType &           dof_handler,
894     const MGLevelObject<VectorType> &data,
895     const std::vector<std::string> & names,
896     const std::vector<DataComponentInterpretation::DataComponentInterpretation>
897       &data_component_interpretation = std::vector<
898         DataComponentInterpretation::DataComponentInterpretation>());
899 
900   /**
901    * Scalar version of the function above.
902    */
903   template <class VectorType>
904   void
905   add_mg_data_vector(const DoFHandlerType &           dof_handler,
906                      const MGLevelObject<VectorType> &data,
907                      const std::string &              name);
908 
909   /**
910    * Release the pointers to the data vectors. This allows output of a new set
911    * of vectors without supplying the DoF handler again. Therefore, the
912    * DataOut object can be used in an algebraic context. Note that besides the
913    * data vectors also the patches already computed are deleted.
914    */
915   void
916   clear_data_vectors();
917 
918   /**
919    * Release pointers to all input data elements, i.e. pointers to data
920    * vectors and to the DoF handler object. This function may be useful when
921    * you have called the @p build_patches function of derived class, since
922    * then the patches are built and the input data is no more needed, nor is
923    * there a need to reference it. You can then output the patches detached
924    * from the main thread and need not make sure anymore that the DoF handler
925    * object and vectors must not be deleted before the output thread is
926    * finished.
927    */
928   void
929   clear_input_data_references();
930 
931   /**
932    * This function can be used to merge the patches that were created using
933    * the @p build_patches function of the object given as argument into the
934    * list of patches created by this object. This is sometimes handy if one
935    * has, for example, a domain decomposition algorithm where each block is
936    * represented by a DoFHandler of its own, but one wants to output the
937    * solution on all the blocks at the same time.
938    *
939    * For this to work, the given argument and this object need to have the
940    * same number of output vectors, and they need to use the same number of
941    * subdivisions per patch. The output will probably look rather funny if
942    * patches in both objects overlap in space.
943    *
944    * If you call build_patches() for this object after merging in patches, the
945    * previous state is overwritten, and the merged-in patches are lost.
946    *
947    * The second parameter allows to shift each node of the patches in the
948    * object passed in in the first parameter by a certain amount. This is
949    * sometimes useful to generate "exploded" views of a collection of blocks.
950    *
951    * This function will fail if either this or the other object did not yet
952    * set up any patches.
953    */
954   template <typename DoFHandlerType2>
955   void
956   merge_patches(
957     const DataOut_DoFData<DoFHandlerType2, patch_dim, patch_space_dim> &source,
958     const Point<patch_space_dim> &shift = Point<patch_space_dim>());
959 
960   /**
961    * Release the pointers to the data vectors and the DoF handler. You have to
962    * set all data entries again using the add_data_vector() function. The
963    * pointer to the dof handler is cleared as well, along with all other data.
964    * In effect, this function resets everything to a virgin state.
965    */
966   virtual void
967   clear();
968 
969   /**
970    * Determine an estimate for the memory consumption (in bytes) of this
971    * object.
972    */
973   std::size_t
974   memory_consumption() const;
975 
976 protected:
977   /**
978    * Abbreviate the somewhat lengthy name for the Patch class.
979    */
980   using Patch = dealii::DataOutBase::Patch<patch_dim, patch_space_dim>;
981 
982   /**
983    * Pointer to the triangulation object.
984    */
985   SmartPointer<const Triangulation<DoFHandlerType::dimension,
986                                    DoFHandlerType::space_dimension>>
987     triangulation;
988 
989   /**
990    * Pointer to the optional handler object.
991    */
992   SmartPointer<const DoFHandlerType> dofs;
993 
994   /**
995    * List of data elements with vectors of values for each degree of freedom.
996    */
997   std::vector<std::shared_ptr<
998     internal::DataOutImplementation::DataEntryBase<DoFHandlerType>>>
999     dof_data;
1000 
1001   /**
1002    * List of data elements with vectors of values for each cell.
1003    */
1004   std::vector<std::shared_ptr<
1005     internal::DataOutImplementation::DataEntryBase<DoFHandlerType>>>
1006     cell_data;
1007 
1008   /**
1009    * This is a list of patches that is created each time build_patches() is
1010    * called. These patches are used in the output routines of the base
1011    * classes.
1012    */
1013   std::vector<Patch> patches;
1014 
1015   /**
1016    * Function by which the base class's functions get to know what patches
1017    * they shall write to a file.
1018    */
1019   virtual const std::vector<Patch> &
1020   get_patches() const override;
1021 
1022   /**
1023    * Virtual function through which the names of data sets are obtained by the
1024    * output functions of the base class.
1025    */
1026   virtual std::vector<std::string>
1027   get_dataset_names() const override;
1028 
1029   /**
1030    * Extracts the finite elements stored in the dof_data object, including a
1031    * dummy object of FE_DGQ<dim>(0) in case only the triangulation is used.
1032    */
1033   std::vector<
1034     std::shared_ptr<dealii::hp::FECollection<DoFHandlerType::dimension,
1035                                              DoFHandlerType::space_dimension>>>
1036   get_fes() const;
1037 
1038   /**
1039    * Overload of the respective DataOutInterface::get_nonscalar_data_ranges()
1040    * function. See there for a more extensive documentation.
1041    */
1042   virtual std::vector<
1043     std::tuple<unsigned int,
1044                unsigned int,
1045                std::string,
1046                DataComponentInterpretation::DataComponentInterpretation>>
1047   get_nonscalar_data_ranges() const override;
1048 
1049   // Make all template siblings friends. Needed for the merge_patches()
1050   // function.
1051   template <class, int, int>
1052   friend class DataOut_DoFData;
1053 
1054   /**
1055    */
1056   template <int, class>
1057   friend class MGDataOut;
1058 
1059 private:
1060   /**
1061    * Common function called by the four public add_data_vector methods.
1062    */
1063   template <class VectorType>
1064   void
1065   add_data_vector_internal(
1066     const DoFHandlerType *          dof_handler,
1067     const VectorType &              data,
1068     const std::vector<std::string> &names,
1069     const DataVectorType            type,
1070     const std::vector<DataComponentInterpretation::DataComponentInterpretation>
1071       &        data_component_interpretation,
1072     const bool deduce_output_names);
1073 };
1074 
1075 
1076 
1077 // -------------------- template and inline functions ------------------------
1078 template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
1079 template <typename VectorType>
1080 void
add_data_vector(const VectorType & vec,const std::string & name,const DataVectorType type,const std::vector<DataComponentInterpretation::DataComponentInterpretation> & data_component_interpretation)1081 DataOut_DoFData<DoFHandlerType, patch_dim, patch_space_dim>::add_data_vector(
1082   const VectorType &   vec,
1083   const std::string &  name,
1084   const DataVectorType type,
1085   const std::vector<DataComponentInterpretation::DataComponentInterpretation>
1086     &data_component_interpretation)
1087 {
1088   Assert(triangulation != nullptr,
1089          Exceptions::DataOutImplementation::ExcNoTriangulationSelected());
1090   std::vector<std::string> names(1, name);
1091   add_data_vector_internal(
1092     dofs, vec, names, type, data_component_interpretation, true);
1093 }
1094 
1095 
1096 
1097 template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
1098 template <typename VectorType>
1099 void
add_data_vector(const VectorType & vec,const std::vector<std::string> & names,const DataVectorType type,const std::vector<DataComponentInterpretation::DataComponentInterpretation> & data_component_interpretation)1100 DataOut_DoFData<DoFHandlerType, patch_dim, patch_space_dim>::add_data_vector(
1101   const VectorType &              vec,
1102   const std::vector<std::string> &names,
1103   const DataVectorType            type,
1104   const std::vector<DataComponentInterpretation::DataComponentInterpretation>
1105     &data_component_interpretation)
1106 {
1107   Assert(triangulation != nullptr,
1108          Exceptions::DataOutImplementation::ExcNoTriangulationSelected());
1109   add_data_vector_internal(
1110     dofs, vec, names, type, data_component_interpretation, false);
1111 }
1112 
1113 
1114 
1115 template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
1116 template <typename VectorType>
1117 void
add_data_vector(const DoFHandlerType & dof_handler,const VectorType & data,const std::string & name,const std::vector<DataComponentInterpretation::DataComponentInterpretation> & data_component_interpretation)1118 DataOut_DoFData<DoFHandlerType, patch_dim, patch_space_dim>::add_data_vector(
1119   const DoFHandlerType &dof_handler,
1120   const VectorType &    data,
1121   const std::string &   name,
1122   const std::vector<DataComponentInterpretation::DataComponentInterpretation>
1123     &data_component_interpretation)
1124 {
1125   std::vector<std::string> names(1, name);
1126   add_data_vector_internal(&dof_handler,
1127                            data,
1128                            names,
1129                            type_dof_data,
1130                            data_component_interpretation,
1131                            true);
1132 }
1133 
1134 
1135 
1136 template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
1137 template <typename VectorType>
1138 void
add_data_vector(const DoFHandlerType & dof_handler,const VectorType & data,const std::vector<std::string> & names,const std::vector<DataComponentInterpretation::DataComponentInterpretation> & data_component_interpretation)1139 DataOut_DoFData<DoFHandlerType, patch_dim, patch_space_dim>::add_data_vector(
1140   const DoFHandlerType &          dof_handler,
1141   const VectorType &              data,
1142   const std::vector<std::string> &names,
1143   const std::vector<DataComponentInterpretation::DataComponentInterpretation>
1144     &data_component_interpretation)
1145 {
1146   add_data_vector_internal(&dof_handler,
1147                            data,
1148                            names,
1149                            type_dof_data,
1150                            data_component_interpretation,
1151                            false);
1152 }
1153 
1154 
1155 
1156 template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
1157 template <typename VectorType>
1158 void
add_data_vector(const VectorType & vec,const DataPostprocessor<DoFHandlerType::space_dimension> & data_postprocessor)1159 DataOut_DoFData<DoFHandlerType, patch_dim, patch_space_dim>::add_data_vector(
1160   const VectorType &                                        vec,
1161   const DataPostprocessor<DoFHandlerType::space_dimension> &data_postprocessor)
1162 {
1163   Assert(dofs != nullptr,
1164          Exceptions::DataOutImplementation::ExcNoDoFHandlerSelected());
1165   add_data_vector(*dofs, vec, data_postprocessor);
1166 }
1167 
1168 
1169 
1170 template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
1171 template <typename DoFHandlerType2>
1172 void
merge_patches(const DataOut_DoFData<DoFHandlerType2,patch_dim,patch_space_dim> & source,const Point<patch_space_dim> & shift)1173 DataOut_DoFData<DoFHandlerType, patch_dim, patch_space_dim>::merge_patches(
1174   const DataOut_DoFData<DoFHandlerType2, patch_dim, patch_space_dim> &source,
1175   const Point<patch_space_dim> &                                      shift)
1176 {
1177   const std::vector<Patch> &source_patches = source.get_patches();
1178   Assert((patches.size() != 0) && (source_patches.size() != 0),
1179          ExcMessage("When calling this function, both the current "
1180                     "object and the one being merged need to have a "
1181                     "nonzero number of patches associated with it. "
1182                     "Either you called this function on objects that "
1183                     "are empty, or you may have forgotten to call "
1184                     "the 'build_patches()' function."));
1185   // Check equality of component names
1186   Assert(get_dataset_names() == source.get_dataset_names(),
1187          Exceptions::DataOutImplementation::ExcIncompatibleDatasetNames());
1188 
1189   // Make sure patches are compatible. Ideally, we would check that all input
1190   // patches from both collections are all compatible, but we'll be content
1191   // with checking that just the first ones from both sources are.
1192   //
1193   // We check compatibility by testing that both sets of patches result
1194   // from the same number of subdivisions, and that they have the same
1195   // number of source vectors (they really should, since we already checked
1196   // that there are the same number of source components above, but you
1197   // never know). This implies that the data should have the same number of
1198   // columns. They should really have the same number of rows as well,
1199   // but depending on whether a patch has points included or not, the
1200   // number of rows may or may not include coordinates for the points,
1201   // and the comparison has to account for that because in each source
1202   // stream, the patches may include some that have points included.
1203   Assert(patches[0].n_subdivisions == source_patches[0].n_subdivisions,
1204          Exceptions::DataOutImplementation::ExcIncompatiblePatchLists());
1205   Assert(patches[0].data.n_cols() == source_patches[0].data.n_cols(),
1206          Exceptions::DataOutImplementation::ExcIncompatiblePatchLists());
1207   Assert((patches[0].data.n_rows() +
1208           (patches[0].points_are_available ? 0 : patch_space_dim)) ==
1209            (source_patches[0].data.n_rows() +
1210             (source_patches[0].points_are_available ? 0 : patch_space_dim)),
1211          Exceptions::DataOutImplementation::ExcIncompatiblePatchLists());
1212 
1213   // check equality of the vector data
1214   // specifications
1215   Assert(get_nonscalar_data_ranges().size() ==
1216            source.get_nonscalar_data_ranges().size(),
1217          ExcMessage("Both sources need to declare the same components "
1218                     "as vectors."));
1219   for (unsigned int i = 0; i < get_nonscalar_data_ranges().size(); ++i)
1220     {
1221       Assert(std::get<0>(get_nonscalar_data_ranges()[i]) ==
1222                std::get<0>(source.get_nonscalar_data_ranges()[i]),
1223              ExcMessage("Both sources need to declare the same components "
1224                         "as vectors."));
1225       Assert(std::get<1>(get_nonscalar_data_ranges()[i]) ==
1226                std::get<1>(source.get_nonscalar_data_ranges()[i]),
1227              ExcMessage("Both sources need to declare the same components "
1228                         "as vectors."));
1229       Assert(std::get<2>(get_nonscalar_data_ranges()[i]) ==
1230                std::get<2>(source.get_nonscalar_data_ranges()[i]),
1231              ExcMessage("Both sources need to declare the same components "
1232                         "as vectors."));
1233     }
1234 
1235   // merge patches. store old number
1236   // of elements, since we need to
1237   // adjust patch numbers, etc
1238   // afterwards
1239   const unsigned int old_n_patches = patches.size();
1240   patches.insert(patches.end(), source_patches.begin(), source_patches.end());
1241 
1242   // perform shift, if so desired
1243   if (shift != Point<patch_space_dim>())
1244     for (unsigned int i = old_n_patches; i < patches.size(); ++i)
1245       for (const unsigned int v : GeometryInfo<patch_dim>::vertex_indices())
1246         patches[i].vertices[v] += shift;
1247 
1248   // adjust patch numbers
1249   for (unsigned int i = old_n_patches; i < patches.size(); ++i)
1250     patches[i].patch_index += old_n_patches;
1251 
1252   // adjust patch neighbors
1253   for (unsigned int i = old_n_patches; i < patches.size(); ++i)
1254     for (const unsigned int n : GeometryInfo<patch_dim>::face_indices())
1255       if (patches[i].neighbors[n] != Patch::no_neighbor)
1256         patches[i].neighbors[n] += old_n_patches;
1257 }
1258 
1259 
1260 DEAL_II_NAMESPACE_CLOSE
1261 
1262 #endif
1263