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