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_templates_h
17 #define dealii_data_out_dof_data_templates_h
18
19
20 #include <deal.II/base/config.h>
21
22 #include <deal.II/base/memory_consumption.h>
23 #include <deal.II/base/numbers.h>
24 #include <deal.II/base/quadrature_lib.h>
25 #include <deal.II/base/signaling_nan.h>
26 #include <deal.II/base/utilities.h>
27 #include <deal.II/base/work_stream.h>
28
29 #include <deal.II/dofs/dof_accessor.h>
30 #include <deal.II/dofs/dof_handler.h>
31
32 #include <deal.II/fe/fe_dgq.h>
33 #include <deal.II/fe/fe_values.h>
34 #include <deal.II/fe/mapping.h>
35
36 #include <deal.II/grid/tria.h>
37 #include <deal.II/grid/tria_iterator.h>
38
39 #include <deal.II/hp/dof_handler.h>
40 #include <deal.II/hp/fe_collection.h>
41 #include <deal.II/hp/fe_values.h>
42 #include <deal.II/hp/q_collection.h>
43
44 #include <deal.II/lac/vector.h>
45
46 #include <deal.II/numerics/data_out.h>
47 #include <deal.II/numerics/data_out_dof_data.h>
48
49 #include <deal.II/simplex/fe_lib.h>
50
51 #include <memory>
52 #include <string>
53 #include <utility>
54 #include <vector>
55
56 DEAL_II_NAMESPACE_OPEN
57
58
59 namespace internal
60 {
61 namespace DataOutImplementation
62 {
63 template <int dim, int spacedim>
ParallelDataBase(const unsigned int n_datasets,const unsigned int n_subdivisions,const std::vector<unsigned int> & n_postprocessor_outputs,const dealii::Mapping<dim,spacedim> & mapping,const std::vector<std::shared_ptr<dealii::hp::FECollection<dim,spacedim>>> & finite_elements,const UpdateFlags update_flags,const bool use_face_values)64 ParallelDataBase<dim, spacedim>::ParallelDataBase(
65 const unsigned int n_datasets,
66 const unsigned int n_subdivisions,
67 const std::vector<unsigned int> & n_postprocessor_outputs,
68 const dealii::Mapping<dim, spacedim> &mapping,
69 const std::vector<
70 std::shared_ptr<dealii::hp::FECollection<dim, spacedim>>>
71 & finite_elements,
72 const UpdateFlags update_flags,
73 const bool use_face_values)
74 : ParallelDataBase<dim, spacedim>(
75 n_datasets,
76 n_subdivisions,
77 n_postprocessor_outputs,
78 dealii::hp::MappingCollection<dim, spacedim>(mapping),
79 finite_elements,
80 update_flags,
81 use_face_values)
82 {}
83
84
85 template <int dim, int spacedim>
ParallelDataBase(const unsigned int n_datasets,const unsigned int n_subdivisions,const std::vector<unsigned int> & n_postprocessor_outputs,const dealii::hp::MappingCollection<dim,spacedim> & mapping,const std::vector<std::shared_ptr<dealii::hp::FECollection<dim,spacedim>>> & finite_elements,const UpdateFlags update_flags,const bool use_face_values)86 ParallelDataBase<dim, spacedim>::ParallelDataBase(
87 const unsigned int n_datasets,
88 const unsigned int n_subdivisions,
89 const std::vector<unsigned int> &n_postprocessor_outputs,
90 const dealii::hp::MappingCollection<dim, spacedim> &mapping,
91 const std::vector<
92 std::shared_ptr<dealii::hp::FECollection<dim, spacedim>>>
93 & finite_elements,
94 const UpdateFlags update_flags,
95 const bool use_face_values)
96 : n_datasets(n_datasets)
97 , n_subdivisions(n_subdivisions)
98 , postprocessed_values(n_postprocessor_outputs.size())
99 , mapping_collection(mapping)
100 , finite_elements(finite_elements)
101 , update_flags(update_flags)
102 {
103 unsigned int n_q_points = 0;
104 if (use_face_values == false)
105 {
106 const bool do_simplex =
107 std::any_of(finite_elements.begin(),
108 finite_elements.end(),
109 [](const auto &fe) {
110 for (unsigned int i = 0; i < fe->size(); ++i)
111 if ((*fe)[i].reference_cell_type() !=
112 ReferenceCell::get_hypercube(dim))
113 return true;
114 return false;
115 });
116
117 const bool do_hex =
118 std::any_of(finite_elements.begin(),
119 finite_elements.end(),
120 [](const auto &fe) {
121 for (unsigned int i = 0; i < fe->size(); ++i)
122 if ((*fe)[i].reference_cell_type() ==
123 ReferenceCell::get_hypercube(dim))
124 return true;
125 return false;
126 });
127
128 std::unique_ptr<dealii::Quadrature<dim>> quadrature1;
129 std::unique_ptr<dealii::Quadrature<dim>> quadrature2;
130
131 if (do_simplex)
132 {
133 Assert(1 <= n_subdivisions && n_subdivisions <= 2,
134 ExcNotImplemented());
135 quadrature1.reset(new Quadrature<dim>(
136 Simplex::FE_P<dim, spacedim>(n_subdivisions == 1 ? 1 : 2)
137 .get_unit_support_points()));
138 }
139
140 if (do_hex)
141 {
142 quadrature2.reset(new Quadrature<dim>(
143 QIterated<dim>(QTrapezoid<1>(), n_subdivisions)));
144 }
145
146 n_q_points = std::max(do_simplex ? quadrature1->size() : 0,
147 do_hex ? quadrature2->size() : 1);
148
149 x_fe_values.resize(this->finite_elements.size());
150 for (unsigned int i = 0; i < this->finite_elements.size(); ++i)
151 {
152 // check if there is a finite element that is equal to the present
153 // one, then we can re-use the FEValues object
154 for (unsigned int j = 0; j < i; ++j)
155 if (this->finite_elements[i].get() ==
156 this->finite_elements[j].get())
157 {
158 x_fe_values[i] = x_fe_values[j];
159 break;
160 }
161 if (x_fe_values[i].get() == nullptr)
162 {
163 dealii::hp::QCollection<dim> quadrature;
164
165 for (unsigned int j = 0; j < this->finite_elements[i]->size();
166 ++j)
167 if ((*this->finite_elements[i])[j].reference_cell_type() !=
168 ReferenceCell::get_hypercube(dim))
169 quadrature.push_back(*quadrature1);
170 else
171 quadrature.push_back(*quadrature2);
172
173 x_fe_values[i] =
174 std::make_shared<dealii::hp::FEValues<dim, spacedim>>(
175 this->mapping_collection,
176 *this->finite_elements[i],
177 quadrature,
178 do_simplex ? (update_flags | update_quadrature_points) :
179 update_flags);
180 }
181 }
182 }
183 else
184 {
185 dealii::hp::QCollection<dim - 1> quadrature(
186 QIterated<dim - 1>(QTrapezoid<1>(), n_subdivisions));
187 n_q_points = quadrature[0].size();
188 x_fe_face_values.resize(this->finite_elements.size());
189 for (unsigned int i = 0; i < this->finite_elements.size(); ++i)
190 {
191 // check if there is a finite element that is equal to the present
192 // one, then we can re-use the FEValues object
193 for (unsigned int j = 0; j < i; ++j)
194 if (this->finite_elements[i].get() ==
195 this->finite_elements[j].get())
196 {
197 x_fe_face_values[i] = x_fe_face_values[j];
198 break;
199 }
200 if (x_fe_face_values[i].get() == nullptr)
201 x_fe_face_values[i] =
202 std::make_shared<dealii::hp::FEFaceValues<dim, spacedim>>(
203 this->mapping_collection,
204 *this->finite_elements[i],
205 quadrature,
206 this->update_flags);
207 }
208 }
209
210 patch_values_scalar.solution_values.resize(n_q_points);
211 patch_values_scalar.solution_gradients.resize(n_q_points);
212 patch_values_scalar.solution_hessians.resize(n_q_points);
213 patch_values_system.solution_values.resize(n_q_points);
214 patch_values_system.solution_gradients.resize(n_q_points);
215 patch_values_system.solution_hessians.resize(n_q_points);
216
217 for (unsigned int dataset = 0; dataset < n_postprocessor_outputs.size();
218 ++dataset)
219 if (n_postprocessor_outputs[dataset] != 0)
220 postprocessed_values[dataset].resize(
221 n_q_points,
222 dealii::Vector<double>(n_postprocessor_outputs[dataset]));
223 }
224
225
226
227 // implement copy constructor to create a thread's own version of
228 // x_fe_values
229 template <int dim, int spacedim>
ParallelDataBase(const ParallelDataBase<dim,spacedim> & data)230 ParallelDataBase<dim, spacedim>::ParallelDataBase(
231 const ParallelDataBase<dim, spacedim> &data)
232 : n_datasets(data.n_datasets)
233 , n_subdivisions(data.n_subdivisions)
234 , patch_values_scalar(data.patch_values_scalar)
235 , patch_values_system(data.patch_values_system)
236 , postprocessed_values(data.postprocessed_values)
237 , mapping_collection(data.mapping_collection)
238 , finite_elements(data.finite_elements)
239 , update_flags(data.update_flags)
240 {
241 if (data.x_fe_values.empty() == false)
242 {
243 Assert(data.x_fe_face_values.empty() == true, ExcInternalError());
244
245 const bool do_simplex =
246 std::any_of(finite_elements.begin(),
247 finite_elements.end(),
248 [](const auto &fe) {
249 for (unsigned int i = 0; i < fe->size(); ++i)
250 if ((*fe)[i].reference_cell_type() !=
251 ReferenceCell::get_hypercube(dim))
252 return true;
253 return false;
254 });
255
256 const bool do_hex =
257 std::any_of(finite_elements.begin(),
258 finite_elements.end(),
259 [](const auto &fe) {
260 for (unsigned int i = 0; i < fe->size(); ++i)
261 if ((*fe)[i].reference_cell_type() ==
262 ReferenceCell::get_hypercube(dim))
263 return true;
264 return false;
265 });
266
267 std::unique_ptr<dealii::Quadrature<dim>> quadrature1;
268 std::unique_ptr<dealii::Quadrature<dim>> quadrature2;
269
270 if (do_simplex)
271 {
272 Assert(1 <= n_subdivisions && n_subdivisions <= 2,
273 ExcNotImplemented());
274 quadrature1.reset(new Quadrature<dim>(
275 Simplex::FE_P<dim, spacedim>(n_subdivisions == 1 ? 1 : 2)
276 .get_unit_support_points()));
277 }
278
279 if (do_hex)
280 {
281 quadrature2.reset(new Quadrature<dim>(
282 QIterated<dim>(QTrapezoid<1>(), n_subdivisions)));
283 }
284
285 x_fe_values.resize(this->finite_elements.size());
286
287 for (unsigned int i = 0; i < this->finite_elements.size(); ++i)
288 {
289 // check if there is a finite element that is equal to the present
290 // one, then we can re-use the FEValues object
291 for (unsigned int j = 0; j < i; ++j)
292 if (this->finite_elements[i].get() ==
293 this->finite_elements[j].get())
294 {
295 x_fe_values[i] = x_fe_values[j];
296 break;
297 }
298 if (x_fe_values[i].get() == nullptr)
299 {
300 dealii::hp::QCollection<dim> quadrature;
301
302 for (unsigned int j = 0; j < this->finite_elements[i]->size();
303 ++j)
304 if ((*this->finite_elements[i])[j].reference_cell_type() !=
305 ReferenceCell::get_hypercube(dim))
306 quadrature.push_back(*quadrature1);
307 else
308 quadrature.push_back(*quadrature2);
309
310 x_fe_values[i] =
311 std::make_shared<dealii::hp::FEValues<dim, spacedim>>(
312 this->mapping_collection,
313 *this->finite_elements[i],
314 quadrature,
315 do_simplex ? (update_flags | update_quadrature_points) :
316 update_flags);
317 }
318 }
319 }
320 else
321 {
322 dealii::hp::QCollection<dim - 1> quadrature(
323 QIterated<dim - 1>(QTrapezoid<1>(), n_subdivisions));
324 x_fe_face_values.resize(this->finite_elements.size());
325 for (unsigned int i = 0; i < this->finite_elements.size(); ++i)
326 {
327 // check if there is a finite element that is equal to the present
328 // one, then we can re-use the FEValues object
329 for (unsigned int j = 0; j < i; ++j)
330 if (this->finite_elements[i].get() ==
331 this->finite_elements[j].get())
332 {
333 x_fe_face_values[i] = x_fe_face_values[j];
334 break;
335 }
336 if (x_fe_face_values[i].get() == nullptr)
337 x_fe_face_values[i] =
338 std::make_shared<dealii::hp::FEFaceValues<dim, spacedim>>(
339 this->mapping_collection,
340 *this->finite_elements[i],
341 quadrature,
342 this->update_flags);
343 }
344 }
345 }
346
347
348
349 template <int dim, int spacedim>
350 template <typename DoFHandlerType>
351 void
reinit_all_fe_values(std::vector<std::shared_ptr<DataEntryBase<DoFHandlerType>>> & dof_data,const typename dealii::Triangulation<dim,spacedim>::cell_iterator & cell,const unsigned int face)352 ParallelDataBase<dim, spacedim>::reinit_all_fe_values(
353 std::vector<std::shared_ptr<DataEntryBase<DoFHandlerType>>> &dof_data,
354 const typename dealii::Triangulation<dim, spacedim>::cell_iterator &cell,
355 const unsigned int face)
356 {
357 for (unsigned int dataset = 0; dataset < dof_data.size(); ++dataset)
358 {
359 const bool is_duplicate = std::any_of(
360 finite_elements.cbegin(),
361 finite_elements.cbegin() + dataset,
362 [&](const std::shared_ptr<dealii::hp::FECollection<dim, spacedim>>
363 &fe) { return finite_elements[dataset].get() == fe.get(); });
364 if (is_duplicate == false)
365 {
366 if (cell->active())
367 {
368 typename DoFHandlerType::active_cell_iterator dh_cell(
369 &cell->get_triangulation(),
370 cell->level(),
371 cell->index(),
372 dof_data[dataset]->dof_handler);
373 if (x_fe_values.empty())
374 {
375 AssertIndexRange(face, GeometryInfo<dim>::faces_per_cell);
376 x_fe_face_values[dataset]->reinit(dh_cell, face);
377 }
378 else
379 x_fe_values[dataset]->reinit(dh_cell);
380 }
381 else
382 x_fe_values[dataset]->reinit(cell);
383 }
384 }
385 if (dof_data.empty())
386 {
387 if (x_fe_values.empty())
388 {
389 AssertIndexRange(face, GeometryInfo<dim>::faces_per_cell);
390 x_fe_face_values[0]->reinit(cell, face);
391 }
392 else
393 x_fe_values[0]->reinit(cell);
394 }
395 }
396
397
398
399 template <int dim, int spacedim>
400 const FEValuesBase<dim, spacedim> &
get_present_fe_values(const unsigned int dataset)401 ParallelDataBase<dim, spacedim>::get_present_fe_values(
402 const unsigned int dataset) const
403 {
404 AssertIndexRange(dataset, finite_elements.size());
405 if (x_fe_values.empty())
406 return x_fe_face_values[dataset]->get_present_fe_values();
407 else
408 return x_fe_values[dataset]->get_present_fe_values();
409 }
410
411
412
413 template <int dim, int spacedim>
414 void
resize_system_vectors(const unsigned int n_components)415 ParallelDataBase<dim, spacedim>::resize_system_vectors(
416 const unsigned int n_components)
417 {
418 Assert(patch_values_system.solution_values.size() > 0,
419 ExcInternalError());
420 AssertDimension(patch_values_system.solution_values.size(),
421 patch_values_system.solution_gradients.size());
422 AssertDimension(patch_values_system.solution_values.size(),
423 patch_values_system.solution_hessians.size());
424 if (patch_values_system.solution_values[0].size() == n_components)
425 return;
426 for (unsigned int k = 0; k < patch_values_system.solution_values.size();
427 ++k)
428 {
429 patch_values_system.solution_values[k].reinit(n_components);
430 patch_values_system.solution_gradients[k].resize(n_components);
431 patch_values_system.solution_hessians[k].resize(n_components);
432 }
433 }
434
435
436
437 /**
438 * In a WorkStream context, use this function to append the patch computed
439 * by the parallel stage to the array of patches.
440 */
441 template <int dim, int spacedim>
442 void
append_patch_to_list(const DataOutBase::Patch<dim,spacedim> & patch,std::vector<DataOutBase::Patch<dim,spacedim>> & patches)443 append_patch_to_list(
444 const DataOutBase::Patch<dim, spacedim> & patch,
445 std::vector<DataOutBase::Patch<dim, spacedim>> &patches)
446 {
447 patches.push_back(patch);
448 patches.back().patch_index = patches.size() - 1;
449 }
450 } // namespace DataOutImplementation
451 } // namespace internal
452
453 namespace internal
454 {
455 namespace DataOutImplementation
456 {
457 /**
458 * Extract the specified component of a number. This template is used when
459 * the given value is assumed to be a real scalar, so asking for the real
460 * part is the only valid choice for the second argument.
461 */
462 template <typename NumberType>
463 double
get_component(const NumberType value,const ComponentExtractor extract_component)464 get_component(const NumberType value,
465 const ComponentExtractor extract_component)
466 {
467 (void)extract_component;
468 static_assert(
469 numbers::NumberTraits<NumberType>::is_complex == false,
470 "This function must not be called for complex-valued data types.");
471 Assert(extract_component == ComponentExtractor::real_part,
472 ExcMessage("You cannot extract anything other than the real "
473 "part from a real number."));
474 return value;
475 }
476
477
478
479 /**
480 * Extract the specified component of a number. This template is used when
481 * the given value is a complex number
482 */
483 template <typename NumberType>
484 double
get_component(const std::complex<NumberType> & value,const ComponentExtractor extract_component)485 get_component(const std::complex<NumberType> &value,
486 const ComponentExtractor extract_component)
487 {
488 switch (extract_component)
489 {
490 case ComponentExtractor::real_part:
491 return value.real();
492
493 case ComponentExtractor::imaginary_part:
494 return value.imag();
495
496 default:
497 Assert(false, ExcInternalError());
498 }
499
500 return numbers::signaling_nan<double>();
501 }
502
503
504
505 template <int rank, int dim, typename NumberType>
506 Tensor<rank, dim>
get_component(const Tensor<rank,dim,NumberType> & value,const ComponentExtractor extract_component)507 get_component(const Tensor<rank, dim, NumberType> &value,
508 const ComponentExtractor extract_component)
509 {
510 Assert(extract_component == ComponentExtractor::real_part,
511 ExcMessage("You cannot extract anything other than the real "
512 "part from a real number."));
513
514 Tensor<rank, dim, double> t;
515 for (unsigned int d = 0; d < dim; ++d)
516 t[d] = get_component(value[d], extract_component);
517
518 return t;
519 }
520
521
522 /**
523 * Helper class templated on vector type to allow different implementations
524 * to extract information from a vector.
525 */
526 template <typename VectorType>
527 struct VectorHelper
528 {
529 /**
530 * extract the @p indices from @p vector and put them into @p values.
531 */
532 static void
533 extract(const VectorType & vector,
534 const std::vector<types::global_dof_index> &indices,
535 const ComponentExtractor extract_component,
536 std::vector<double> & values);
537 };
538
539
540
541 template <typename VectorType>
542 void
extract(const VectorType & vector,const std::vector<types::global_dof_index> & indices,const ComponentExtractor extract_component,std::vector<double> & values)543 VectorHelper<VectorType>::extract(
544 const VectorType & vector,
545 const std::vector<types::global_dof_index> &indices,
546 const ComponentExtractor extract_component,
547 std::vector<double> & values)
548 {
549 for (unsigned int i = 0; i < values.size(); ++i)
550 values[i] = get_component(vector[indices[i]], extract_component);
551 }
552
553
554
555 #ifdef DEAL_II_WITH_TRILINOS
556 template <>
557 inline void
extract(const LinearAlgebra::EpetraWrappers::Vector &,const std::vector<types::global_dof_index> &,const ComponentExtractor,std::vector<double> &)558 VectorHelper<LinearAlgebra::EpetraWrappers::Vector>::extract(
559 const LinearAlgebra::EpetraWrappers::Vector & /*vector*/,
560 const std::vector<types::global_dof_index> & /*indices*/,
561 const ComponentExtractor /*extract_component*/,
562 std::vector<double> & /*values*/)
563 {
564 // TODO: we don't have element access
565 Assert(false, ExcNotImplemented());
566 }
567 #endif
568
569
570
571 #if defined(DEAL_II_TRILINOS_WITH_TPETRA) && defined(DEAL_II_WITH_MPI)
572 template <>
573 inline void
extract(const LinearAlgebra::TpetraWrappers::Vector<double> &,const std::vector<types::global_dof_index> &,const ComponentExtractor,std::vector<double> &)574 VectorHelper<LinearAlgebra::TpetraWrappers::Vector<double>>::extract(
575 const LinearAlgebra::TpetraWrappers::Vector<double> & /*vector*/,
576 const std::vector<types::global_dof_index> & /*indices*/,
577 const ComponentExtractor /*extract_component*/,
578 std::vector<double> & /*values*/)
579 {
580 // TODO: we don't have element access
581 Assert(false, ExcNotImplemented());
582 }
583
584 template <>
585 inline void
extract(const LinearAlgebra::TpetraWrappers::Vector<float> &,const std::vector<types::global_dof_index> &,const ComponentExtractor,std::vector<double> &)586 VectorHelper<LinearAlgebra::TpetraWrappers::Vector<float>>::extract(
587 const LinearAlgebra::TpetraWrappers::Vector<float> & /*vector*/,
588 const std::vector<types::global_dof_index> & /*indices*/,
589 const ComponentExtractor /*extract_component*/,
590 std::vector<double> & /*values*/)
591 {
592 // TODO: we don't have element access
593 Assert(false, ExcNotImplemented());
594 }
595 #endif
596
597
598
599 template <typename DoFHandlerType>
DataEntryBase(const DoFHandlerType * dofs,const std::vector<std::string> & names_in,const std::vector<DataComponentInterpretation::DataComponentInterpretation> & data_component_interpretation)600 DataEntryBase<DoFHandlerType>::DataEntryBase(
601 const DoFHandlerType * dofs,
602 const std::vector<std::string> &names_in,
603 const std::vector<
604 DataComponentInterpretation::DataComponentInterpretation>
605 &data_component_interpretation)
606 : dof_handler(dofs,
607 typeid(
608 dealii::DataOut_DoFData<DoFHandlerType,
609 DoFHandlerType::dimension,
610 DoFHandlerType::space_dimension>)
611 .name())
612 , names(names_in)
613 , data_component_interpretation(data_component_interpretation)
614 , postprocessor(nullptr, typeid(*this).name())
615 , n_output_variables(names.size())
616 {
617 Assert(names.size() == data_component_interpretation.size(),
618 ExcDimensionMismatch(data_component_interpretation.size(),
619 names.size()));
620
621 // check that the names use only allowed characters
622 for (const auto &name : names)
623 {
624 (void)name;
625 Assert(name.find_first_not_of("abcdefghijklmnopqrstuvwxyz"
626 "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
627 "0123456789_<>()") == std::string::npos,
628 Exceptions::DataOutImplementation::ExcInvalidCharacter(
629 name,
630 name.find_first_not_of("abcdefghijklmnopqrstuvwxyz"
631 "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
632 "0123456789_<>()")));
633 }
634 }
635
636
637
638 template <typename DoFHandlerType>
DataEntryBase(const DoFHandlerType * dofs,const DataPostprocessor<DoFHandlerType::space_dimension> * data_postprocessor)639 DataEntryBase<DoFHandlerType>::DataEntryBase(
640 const DoFHandlerType *dofs,
641 const DataPostprocessor<DoFHandlerType::space_dimension>
642 *data_postprocessor)
643 : dof_handler(dofs,
644 typeid(
645 dealii::DataOut_DoFData<DoFHandlerType,
646 DoFHandlerType::dimension,
647 DoFHandlerType::space_dimension>)
648 .name())
649 , names(data_postprocessor->get_names())
650 , data_component_interpretation(
651 data_postprocessor->get_data_component_interpretation())
652 , postprocessor(data_postprocessor, typeid(*this).name())
653 , n_output_variables(names.size())
654 {
655 Assert(data_postprocessor->get_names().size() ==
656 data_postprocessor->get_data_component_interpretation().size(),
657 ExcDimensionMismatch(
658 data_postprocessor->get_names().size(),
659 data_postprocessor->get_data_component_interpretation().size()));
660
661 // check that the names use only allowed characters
662 for (const auto &name : names)
663 {
664 (void)name;
665 Assert(name.find_first_not_of("abcdefghijklmnopqrstuvwxyz"
666 "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
667 "0123456789_<>()") == std::string::npos,
668 Exceptions::DataOutImplementation::ExcInvalidCharacter(
669 name,
670 name.find_first_not_of("abcdefghijklmnopqrstuvwxyz"
671 "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
672 "0123456789_<>()")));
673 }
674 }
675
676
677
678 /**
679 * Class that stores a pointer to a vector of type equal to the template
680 * argument, and provides the functions to extract data from it.
681 */
682 template <typename DoFHandlerType, typename VectorType>
683 class DataEntry : public DataEntryBase<DoFHandlerType>
684 {
685 public:
686 /**
687 * Constructor. Give a list of names for the individual components of
688 * the vector and their interpretation as scalar or vector data. This
689 * constructor assumes that no postprocessor is going to be used.
690 */
691 DataEntry(const DoFHandlerType * dofs,
692 const VectorType * data,
693 const std::vector<std::string> &names,
694 const std::vector<
695 DataComponentInterpretation::DataComponentInterpretation>
696 &data_component_interpretation);
697
698 /**
699 * Constructor when a data postprocessor is going to be used. In that
700 * case, the names and vector declarations are going to be acquired from
701 * the postprocessor.
702 */
703 DataEntry(const DoFHandlerType *dofs,
704 const VectorType * data,
705 const DataPostprocessor<DoFHandlerType::space_dimension>
706 *data_postprocessor);
707
708 /**
709 * Assuming that the stored vector is a cell vector, extract the given
710 * element from it.
711 */
712 virtual double
713 get_cell_data_value(
714 const unsigned int cell_number,
715 const ComponentExtractor extract_component) const override;
716
717 /**
718 * Given a FEValuesBase object, extract the values on the present cell
719 * from the vector we actually store.
720 */
721 virtual void
722 get_function_values(
723 const FEValuesBase<DoFHandlerType::dimension,
724 DoFHandlerType::space_dimension> &fe_patch_values,
725 const ComponentExtractor extract_component,
726 std::vector<double> &patch_values) const override;
727
728 /**
729 * Given a FEValuesBase object, extract the values on the present cell
730 * from the vector we actually store. This function does the same as the
731 * one above but for vector-valued finite elements.
732 */
733 virtual void
734 get_function_values(
735 const FEValuesBase<DoFHandlerType::dimension,
736 DoFHandlerType::space_dimension> &fe_patch_values,
737 const ComponentExtractor extract_component,
738 std::vector<dealii::Vector<double>> &patch_values_system)
739 const override;
740
741 /**
742 * Given a FEValuesBase object, extract the gradients on the present
743 * cell from the vector we actually store.
744 */
745 virtual void
746 get_function_gradients(
747 const FEValuesBase<DoFHandlerType::dimension,
748 DoFHandlerType::space_dimension> &fe_patch_values,
749 const ComponentExtractor extract_component,
750 std::vector<Tensor<1, DoFHandlerType::space_dimension>>
751 &patch_gradients) const override;
752
753 /**
754 * Given a FEValuesBase object, extract the gradients on the present
755 * cell from the vector we actually store. This function does the same
756 * as the one above but for vector-valued finite elements.
757 */
758 virtual void
759 get_function_gradients(
760 const FEValuesBase<DoFHandlerType::dimension,
761 DoFHandlerType::space_dimension> &fe_patch_values,
762 const ComponentExtractor extract_component,
763 std::vector<std::vector<Tensor<1, DoFHandlerType::space_dimension>>>
764 &patch_gradients_system) const override;
765
766 /**
767 * Given a FEValuesBase object, extract the second derivatives on the
768 * present cell from the vector we actually store.
769 */
770 virtual void
771 get_function_hessians(
772 const FEValuesBase<DoFHandlerType::dimension,
773 DoFHandlerType::space_dimension> &fe_patch_values,
774 const ComponentExtractor extract_component,
775 std::vector<Tensor<2, DoFHandlerType::space_dimension>> &patch_hessians)
776 const override;
777
778 /**
779 * Given a FEValuesBase object, extract the second derivatives on the
780 * present cell from the vector we actually store. This function does
781 * the same as the one above but for vector-valued finite elements.
782 */
783 virtual void
784 get_function_hessians(
785 const FEValuesBase<DoFHandlerType::dimension,
786 DoFHandlerType::space_dimension> &fe_patch_values,
787 const ComponentExtractor extract_component,
788 std::vector<std::vector<Tensor<2, DoFHandlerType::space_dimension>>>
789 &patch_hessians_system) const override;
790
791 /**
792 * Return whether the data represented by (a derived class of) this object
793 * represents a complex-valued (as opposed to real-valued) information.
794 */
795 virtual bool
796 is_complex_valued() const override;
797
798 /**
799 * Clear all references to the vectors.
800 */
801 virtual void
802 clear() override;
803
804 /**
805 * Determine an estimate for the memory consumption (in bytes) of this
806 * object.
807 */
808 virtual std::size_t
809 memory_consumption() const override;
810
811 private:
812 /**
813 * Pointer to the data vector. Note that ownership of the vector pointed
814 * to remains with the caller of this class.
815 */
816 const VectorType *vector;
817 };
818
819
820
821 template <typename DoFHandlerType, typename VectorType>
DataEntry(const DoFHandlerType * dofs,const VectorType * data,const std::vector<std::string> & names,const std::vector<DataComponentInterpretation::DataComponentInterpretation> & data_component_interpretation)822 DataEntry<DoFHandlerType, VectorType>::DataEntry(
823 const DoFHandlerType * dofs,
824 const VectorType * data,
825 const std::vector<std::string> &names,
826 const std::vector<
827 DataComponentInterpretation::DataComponentInterpretation>
828 &data_component_interpretation)
829 : DataEntryBase<DoFHandlerType>(dofs,
830 names,
831 data_component_interpretation)
832 , vector(data)
833 {}
834
835
836
837 template <typename DoFHandlerType, typename VectorType>
DataEntry(const DoFHandlerType * dofs,const VectorType * data,const DataPostprocessor<DoFHandlerType::space_dimension> * data_postprocessor)838 DataEntry<DoFHandlerType, VectorType>::DataEntry(
839 const DoFHandlerType *dofs,
840 const VectorType * data,
841 const DataPostprocessor<DoFHandlerType::space_dimension>
842 *data_postprocessor)
843 : DataEntryBase<DoFHandlerType>(dofs, data_postprocessor)
844 , vector(data)
845 {}
846
847
848
849 template <typename DoFHandlerType, typename VectorType>
850 double
get_cell_data_value(const unsigned int cell_number,const ComponentExtractor extract_component)851 DataEntry<DoFHandlerType, VectorType>::get_cell_data_value(
852 const unsigned int cell_number,
853 const ComponentExtractor extract_component) const
854 {
855 return get_component(
856 internal::ElementAccess<VectorType>::get(*vector, cell_number),
857 extract_component);
858 }
859
860
861
862 template <typename DoFHandlerType, typename VectorType>
863 void
get_function_values(const FEValuesBase<DoFHandlerType::dimension,DoFHandlerType::space_dimension> & fe_patch_values,const ComponentExtractor extract_component,std::vector<dealii::Vector<double>> & patch_values_system)864 DataEntry<DoFHandlerType, VectorType>::get_function_values(
865 const FEValuesBase<DoFHandlerType::dimension,
866 DoFHandlerType::space_dimension> &fe_patch_values,
867 const ComponentExtractor extract_component,
868 std::vector<dealii::Vector<double>> &patch_values_system) const
869 {
870 if (typeid(typename VectorType::value_type) == typeid(double))
871 {
872 Assert(extract_component == ComponentExtractor::real_part,
873 ExcMessage("You cannot extract anything other than the real "
874 "part from a real number."));
875
876 fe_patch_values.get_function_values(
877 *vector,
878 // reinterpret output argument type; because of the 'if' statement
879 // above, this is the identity cast whenever the code is executed,
880 // but the cast is necessary to allow compilation even if we don't
881 // get here
882 reinterpret_cast<
883 std::vector<dealii::Vector<typename VectorType::value_type>> &>(
884 patch_values_system));
885 }
886 else
887 {
888 // OK, so we know that the data type stored by the user-provided
889 // vector is not simply a double. In that case, we need to ask the
890 // FEValuesBase object to extract function values for us from the
891 // evaluation points in the provided data type, and then copy them by
892 // hand into the output location.
893 const unsigned int n_components =
894 fe_patch_values.get_fe().n_components();
895 const unsigned int n_eval_points =
896 fe_patch_values.n_quadrature_points;
897
898 std::vector<dealii::Vector<typename VectorType::value_type>> tmp(
899 n_eval_points);
900 for (unsigned int i = 0; i < n_eval_points; i++)
901 tmp[i].reinit(n_components);
902
903 fe_patch_values.get_function_values(*vector, tmp);
904
905 AssertDimension(patch_values_system.size(), n_eval_points);
906 for (unsigned int i = 0; i < n_eval_points; i++)
907 {
908 AssertDimension(patch_values_system[i].size(), n_components);
909
910 for (unsigned int j = 0; j < n_components; ++j)
911 patch_values_system[i](j) =
912 get_component(tmp[i](j), extract_component);
913 }
914 }
915 }
916
917
918
919 template <typename DoFHandlerType, typename VectorType>
920 void
get_function_values(const FEValuesBase<DoFHandlerType::dimension,DoFHandlerType::space_dimension> & fe_patch_values,const ComponentExtractor extract_component,std::vector<double> & patch_values)921 DataEntry<DoFHandlerType, VectorType>::get_function_values(
922 const FEValuesBase<DoFHandlerType::dimension,
923 DoFHandlerType::space_dimension> &fe_patch_values,
924 const ComponentExtractor extract_component,
925 std::vector<double> & patch_values) const
926 {
927 if (typeid(typename VectorType::value_type) == typeid(double))
928 {
929 Assert(extract_component == ComponentExtractor::real_part,
930 ExcMessage("You cannot extract anything other than the real "
931 "part from a real number."));
932
933 fe_patch_values.get_function_values(
934 *vector,
935 // reinterpret output argument type; because of the 'if' statement
936 // above, this is the identity cast whenever the code is executed,
937 // but the cast is necessary to allow compilation even if we don't
938 // get here
939 reinterpret_cast<std::vector<typename VectorType::value_type> &>(
940 patch_values));
941 }
942 else
943 {
944 std::vector<typename VectorType::value_type> tmp(patch_values.size());
945
946 fe_patch_values.get_function_values(*vector, tmp);
947
948 for (unsigned int i = 0; i < tmp.size(); i++)
949 patch_values[i] = get_component(tmp[i], extract_component);
950 }
951 }
952
953
954
955 template <typename DoFHandlerType, typename VectorType>
956 void
get_function_gradients(const FEValuesBase<DoFHandlerType::dimension,DoFHandlerType::space_dimension> & fe_patch_values,const ComponentExtractor extract_component,std::vector<std::vector<Tensor<1,DoFHandlerType::space_dimension>>> & patch_gradients_system)957 DataEntry<DoFHandlerType, VectorType>::get_function_gradients(
958 const FEValuesBase<DoFHandlerType::dimension,
959 DoFHandlerType::space_dimension> &fe_patch_values,
960 const ComponentExtractor extract_component,
961 std::vector<std::vector<Tensor<1, DoFHandlerType::space_dimension>>>
962 &patch_gradients_system) const
963 {
964 if (typeid(typename VectorType::value_type) == typeid(double))
965 {
966 Assert(extract_component == ComponentExtractor::real_part,
967 ExcMessage("You cannot extract anything other than the real "
968 "part from a real number."));
969
970 fe_patch_values.get_function_gradients(
971 *vector,
972 // reinterpret output argument type; because of the 'if' statement
973 // above, this is the identity cast whenever the code is executed,
974 // but the cast is necessary to allow compilation even if we don't
975 // get here
976 reinterpret_cast<std::vector<
977 std::vector<Tensor<1,
978 DoFHandlerType::space_dimension,
979 typename VectorType::value_type>>> &>(
980 patch_gradients_system));
981 }
982 else
983 {
984 // OK, so we know that the data type stored by the user-provided
985 // vector is not simply a double. In that case, we need to ask the
986 // FEValuesBase object to extract function values for us from the
987 // evaluation points in the provided data type, and then copy them by
988 // hand into the output location.
989 const unsigned int n_components =
990 fe_patch_values.get_fe().n_components();
991 const unsigned int n_eval_points =
992 fe_patch_values.n_quadrature_points;
993
994 std::vector<std::vector<Tensor<1,
995 DoFHandlerType::space_dimension,
996 typename VectorType::value_type>>>
997 tmp(n_eval_points);
998 for (unsigned int i = 0; i < n_eval_points; i++)
999 tmp[i].resize(n_components);
1000
1001 fe_patch_values.get_function_gradients(*vector, tmp);
1002
1003 AssertDimension(patch_gradients_system.size(), n_eval_points);
1004 for (unsigned int i = 0; i < n_eval_points; i++)
1005 {
1006 AssertDimension(patch_gradients_system[i].size(), n_components);
1007
1008 for (unsigned int j = 0; j < n_components; j++)
1009 patch_gradients_system[i][j] =
1010 get_component(tmp[i][j], extract_component);
1011 }
1012 }
1013 }
1014
1015
1016
1017 template <typename DoFHandlerType, typename VectorType>
1018 void
get_function_gradients(const FEValuesBase<DoFHandlerType::dimension,DoFHandlerType::space_dimension> & fe_patch_values,const ComponentExtractor extract_component,std::vector<Tensor<1,DoFHandlerType::space_dimension>> & patch_gradients)1019 DataEntry<DoFHandlerType, VectorType>::get_function_gradients(
1020 const FEValuesBase<DoFHandlerType::dimension,
1021 DoFHandlerType::space_dimension> &fe_patch_values,
1022 const ComponentExtractor extract_component,
1023 std::vector<Tensor<1, DoFHandlerType::space_dimension>> &patch_gradients)
1024 const
1025 {
1026 if (typeid(typename VectorType::value_type) == typeid(double))
1027 {
1028 Assert(extract_component == ComponentExtractor::real_part,
1029 ExcMessage("You cannot extract anything other than the real "
1030 "part from a real number."));
1031
1032 fe_patch_values.get_function_gradients(
1033 *vector,
1034 // reinterpret output argument type; because of the 'if' statement
1035 // above, this is the identity cast whenever the code is executed,
1036 // but the cast is necessary to allow compilation even if we don't
1037 // get here
1038 reinterpret_cast<
1039 std::vector<Tensor<1,
1040 DoFHandlerType::space_dimension,
1041 typename VectorType::value_type>> &>(
1042 patch_gradients));
1043 }
1044 else
1045 {
1046 std::vector<Tensor<1,
1047 DoFHandlerType::space_dimension,
1048 typename VectorType::value_type>>
1049 tmp;
1050 tmp.resize(patch_gradients.size());
1051
1052 fe_patch_values.get_function_gradients(*vector, tmp);
1053
1054 for (unsigned int i = 0; i < tmp.size(); i++)
1055 patch_gradients[i] = get_component(tmp[i], extract_component);
1056 }
1057 }
1058
1059
1060
1061 template <typename DoFHandlerType, typename VectorType>
1062 void
get_function_hessians(const FEValuesBase<DoFHandlerType::dimension,DoFHandlerType::space_dimension> & fe_patch_values,const ComponentExtractor extract_component,std::vector<std::vector<Tensor<2,DoFHandlerType::space_dimension>>> & patch_hessians_system)1063 DataEntry<DoFHandlerType, VectorType>::get_function_hessians(
1064 const FEValuesBase<DoFHandlerType::dimension,
1065 DoFHandlerType::space_dimension> &fe_patch_values,
1066 const ComponentExtractor extract_component,
1067 std::vector<std::vector<Tensor<2, DoFHandlerType::space_dimension>>>
1068 &patch_hessians_system) const
1069 {
1070 if (typeid(typename VectorType::value_type) == typeid(double))
1071 {
1072 Assert(extract_component == ComponentExtractor::real_part,
1073 ExcMessage("You cannot extract anything other than the real "
1074 "part from a real number."));
1075
1076 fe_patch_values.get_function_hessians(
1077 *vector,
1078 // reinterpret output argument type; because of the 'if' statement
1079 // above, this is the identity cast whenever the code is executed,
1080 // but the cast is necessary to allow compilation even if we don't
1081 // get here
1082 reinterpret_cast<std::vector<
1083 std::vector<Tensor<2,
1084 DoFHandlerType::space_dimension,
1085 typename VectorType::value_type>>> &>(
1086 patch_hessians_system));
1087 }
1088 else
1089 {
1090 // OK, so we know that the data type stored by the user-provided
1091 // vector is not simply a double. In that case, we need to ask the
1092 // FEValuesBase object to extract function values for us from the
1093 // evaluation points in the provided data type, and then copy them by
1094 // hand into the output location.
1095 const unsigned int n_components =
1096 fe_patch_values.get_fe().n_components();
1097 const unsigned int n_eval_points =
1098 fe_patch_values.n_quadrature_points;
1099
1100 std::vector<std::vector<Tensor<2,
1101 DoFHandlerType::space_dimension,
1102 typename VectorType::value_type>>>
1103 tmp(n_eval_points);
1104 for (unsigned int i = 0; i < n_eval_points; i++)
1105 tmp[i].resize(n_components);
1106
1107 fe_patch_values.get_function_hessians(*vector, tmp);
1108
1109 AssertDimension(patch_hessians_system.size(), n_eval_points);
1110 for (unsigned int i = 0; i < n_eval_points; i++)
1111 {
1112 AssertDimension(patch_hessians_system[i].size(), n_components);
1113
1114 for (unsigned int j = 0; j < n_components; j++)
1115 patch_hessians_system[i][j] =
1116 get_component(tmp[i][j], extract_component);
1117 }
1118 }
1119 }
1120
1121
1122
1123 template <typename DoFHandlerType, typename VectorType>
1124 void
get_function_hessians(const FEValuesBase<DoFHandlerType::dimension,DoFHandlerType::space_dimension> & fe_patch_values,const ComponentExtractor extract_component,std::vector<Tensor<2,DoFHandlerType::space_dimension>> & patch_hessians)1125 DataEntry<DoFHandlerType, VectorType>::get_function_hessians(
1126 const FEValuesBase<DoFHandlerType::dimension,
1127 DoFHandlerType::space_dimension> &fe_patch_values,
1128 const ComponentExtractor extract_component,
1129 std::vector<Tensor<2, DoFHandlerType::space_dimension>> &patch_hessians)
1130 const
1131 {
1132 if (typeid(typename VectorType::value_type) == typeid(double))
1133 {
1134 Assert(extract_component == ComponentExtractor::real_part,
1135 ExcMessage("You cannot extract anything other than the real "
1136 "part from a real number."));
1137
1138 fe_patch_values.get_function_hessians(
1139 *vector,
1140 // reinterpret output argument type; because of the 'if' statement
1141 // above, this is the identity cast whenever the code is executed,
1142 // but the cast is necessary to allow compilation even if we don't
1143 // get here
1144 reinterpret_cast<
1145 std::vector<Tensor<2,
1146 DoFHandlerType::space_dimension,
1147 typename VectorType::value_type>> &>(
1148 patch_hessians));
1149 }
1150 else
1151 {
1152 std::vector<Tensor<2,
1153 DoFHandlerType::space_dimension,
1154 typename VectorType::value_type>>
1155 tmp(patch_hessians.size());
1156
1157 fe_patch_values.get_function_hessians(*vector, tmp);
1158
1159 for (unsigned int i = 0; i < tmp.size(); i++)
1160 patch_hessians[i] = get_component(tmp[i], extract_component);
1161 }
1162 }
1163
1164
1165
1166 template <typename DoFHandlerType, typename VectorType>
1167 bool
is_complex_valued()1168 DataEntry<DoFHandlerType, VectorType>::is_complex_valued() const
1169 {
1170 return numbers::NumberTraits<typename VectorType::value_type>::is_complex;
1171 }
1172
1173
1174
1175 template <typename DoFHandlerType, typename VectorType>
1176 std::size_t
memory_consumption()1177 DataEntry<DoFHandlerType, VectorType>::memory_consumption() const
1178 {
1179 return (sizeof(vector) +
1180 MemoryConsumption::memory_consumption(this->names));
1181 }
1182
1183
1184
1185 template <typename DoFHandlerType, typename VectorType>
1186 void
clear()1187 DataEntry<DoFHandlerType, VectorType>::clear()
1188 {
1189 vector = nullptr;
1190 this->dof_handler = nullptr;
1191 }
1192
1193
1194
1195 /**
1196 * Like DataEntry, but used to look up data from multigrid computations.
1197 * Data will use level-DoF indices to look up in a
1198 * MGLevelObject<VectorType> given on the specific level instead of
1199 * interpolating data to coarser cells.
1200 */
1201 template <typename DoFHandlerType, typename VectorType>
1202 class MGDataEntry : public DataEntryBase<DoFHandlerType>
1203 {
1204 public:
MGDataEntry(const DoFHandlerType * dofs,const MGLevelObject<VectorType> * vectors,const std::vector<std::string> & names,const std::vector<DataComponentInterpretation::DataComponentInterpretation> & data_component_interpretation)1205 MGDataEntry(const DoFHandlerType * dofs,
1206 const MGLevelObject<VectorType> *vectors,
1207 const std::vector<std::string> & names,
1208 const std::vector<
1209 DataComponentInterpretation::DataComponentInterpretation>
1210 &data_component_interpretation)
1211 : DataEntryBase<DoFHandlerType>(dofs,
1212 names,
1213 data_component_interpretation)
1214 , vectors(vectors)
1215 {}
1216
1217 virtual double
1218 get_cell_data_value(
1219 const unsigned int cell_number,
1220 const ComponentExtractor extract_component) const override;
1221
1222 virtual void
1223 get_function_values(
1224 const FEValuesBase<DoFHandlerType::dimension,
1225 DoFHandlerType::space_dimension> &fe_patch_values,
1226 const ComponentExtractor extract_component,
1227 std::vector<double> &patch_values) const override;
1228
1229 /**
1230 * Given a FEValuesBase object, extract the values on the present cell
1231 * from the vector we actually store. This function does the same as the
1232 * one above but for vector-valued finite elements.
1233 */
1234 virtual void
1235 get_function_values(
1236 const FEValuesBase<DoFHandlerType::dimension,
1237 DoFHandlerType::space_dimension> &fe_patch_values,
1238 const ComponentExtractor extract_component,
1239 std::vector<dealii::Vector<double>> &patch_values_system)
1240 const override;
1241
1242 /**
1243 * Given a FEValuesBase object, extract the gradients on the present
1244 * cell from the vector we actually store.
1245 */
1246 virtual void
get_function_gradients(const FEValuesBase<DoFHandlerType::dimension,DoFHandlerType::space_dimension> &,const ComponentExtractor,std::vector<Tensor<1,DoFHandlerType::space_dimension>> &)1247 get_function_gradients(
1248 const FEValuesBase<DoFHandlerType::dimension,
1249 DoFHandlerType::space_dimension>
1250 & /*fe_patch_values*/,
1251 const ComponentExtractor /*extract_component*/,
1252 std::vector<Tensor<1, DoFHandlerType::space_dimension>>
1253 & /*patch_gradients*/) const override
1254 {
1255 Assert(false, ExcNotImplemented());
1256 }
1257
1258 /**
1259 * Given a FEValuesBase object, extract the gradients on the present
1260 * cell from the vector we actually store. This function does the same
1261 * as the one above but for vector-valued finite elements.
1262 */
1263 virtual void
get_function_gradients(const FEValuesBase<DoFHandlerType::dimension,DoFHandlerType::space_dimension> &,const ComponentExtractor,std::vector<std::vector<Tensor<1,DoFHandlerType::space_dimension>>> &)1264 get_function_gradients(
1265 const FEValuesBase<DoFHandlerType::dimension,
1266 DoFHandlerType::space_dimension>
1267 & /*fe_patch_values*/,
1268 const ComponentExtractor /*extract_component*/,
1269 std::vector<std::vector<Tensor<1, DoFHandlerType::space_dimension>>>
1270 & /*patch_gradients_system*/) const override
1271 {
1272 Assert(false, ExcNotImplemented());
1273 }
1274
1275
1276 /**
1277 * Given a FEValuesBase object, extract the second derivatives on the
1278 * present cell from the vector we actually store.
1279 */
1280 virtual void
get_function_hessians(const FEValuesBase<DoFHandlerType::dimension,DoFHandlerType::space_dimension> &,const ComponentExtractor,std::vector<Tensor<2,DoFHandlerType::space_dimension>> &)1281 get_function_hessians(
1282 const FEValuesBase<DoFHandlerType::dimension,
1283 DoFHandlerType::space_dimension>
1284 & /*fe_patch_values*/,
1285 const ComponentExtractor /*extract_component*/,
1286 std::vector<Tensor<2, DoFHandlerType::space_dimension>>
1287 & /*patch_hessians*/) const override
1288 {
1289 Assert(false, ExcNotImplemented());
1290 }
1291
1292 /**
1293 * Given a FEValuesBase object, extract the second derivatives on the
1294 * present cell from the vector we actually store. This function does
1295 * the same as the one above but for vector-valued finite elements.
1296 */
1297 virtual void
get_function_hessians(const FEValuesBase<DoFHandlerType::dimension,DoFHandlerType::space_dimension> &,const ComponentExtractor,std::vector<std::vector<Tensor<2,DoFHandlerType::space_dimension>>> &)1298 get_function_hessians(
1299 const FEValuesBase<DoFHandlerType::dimension,
1300 DoFHandlerType::space_dimension>
1301 & /*fe_patch_values*/,
1302 const ComponentExtractor /*extract_component*/,
1303 std::vector<std::vector<Tensor<2, DoFHandlerType::space_dimension>>>
1304 & /*patch_hessians_system*/) const override
1305 {
1306 Assert(false, ExcNotImplemented());
1307 }
1308
1309 /**
1310 * Return whether the data represented by (a derived class of) this object
1311 * represents a complex-valued (as opposed to real-valued) information.
1312 */
1313 virtual bool
is_complex_valued()1314 is_complex_valued() const override
1315 {
1316 Assert(
1317 numbers::NumberTraits<typename VectorType::value_type>::is_complex ==
1318 false,
1319 ExcNotImplemented());
1320 return numbers::NumberTraits<
1321 typename VectorType::value_type>::is_complex;
1322 }
1323
1324 /**
1325 * Clear all references to the vectors.
1326 */
1327 virtual void
clear()1328 clear() override
1329 {
1330 vectors = nullptr;
1331 }
1332
1333 /**
1334 * Determine an estimate for the memory consumption (in bytes) of this
1335 * object.
1336 */
1337 virtual std::size_t
memory_consumption()1338 memory_consumption() const override
1339 {
1340 return sizeof(vectors);
1341 }
1342
1343 private:
1344 const MGLevelObject<VectorType> *vectors;
1345 };
1346
1347
1348
1349 template <typename DoFHandlerType, typename VectorType>
1350 double
get_cell_data_value(const unsigned int cell_number,const ComponentExtractor extract_component)1351 MGDataEntry<DoFHandlerType, VectorType>::get_cell_data_value(
1352 const unsigned int cell_number,
1353 const ComponentExtractor extract_component) const
1354 {
1355 Assert(false, ExcNotImplemented());
1356
1357 (void)cell_number;
1358 (void)extract_component;
1359 return 0.0;
1360 }
1361
1362
1363
1364 template <typename DoFHandlerType, typename VectorType>
1365 void
get_function_values(const FEValuesBase<DoFHandlerType::dimension,DoFHandlerType::space_dimension> & fe_patch_values,const ComponentExtractor extract_component,std::vector<double> & patch_values)1366 MGDataEntry<DoFHandlerType, VectorType>::get_function_values(
1367 const FEValuesBase<DoFHandlerType::dimension,
1368 DoFHandlerType::space_dimension> &fe_patch_values,
1369 const ComponentExtractor extract_component,
1370 std::vector<double> & patch_values) const
1371 {
1372 Assert(extract_component == ComponentExtractor::real_part,
1373 ExcNotImplemented());
1374
1375 const typename DoFHandlerType::level_cell_iterator dof_cell(
1376 &fe_patch_values.get_cell()->get_triangulation(),
1377 fe_patch_values.get_cell()->level(),
1378 fe_patch_values.get_cell()->index(),
1379 this->dof_handler);
1380
1381 const VectorType *vector = &((*vectors)[dof_cell->level()]);
1382
1383 const unsigned int dofs_per_cell =
1384 this->dof_handler->get_fe(0).n_dofs_per_cell();
1385
1386 std::vector<types::global_dof_index> dof_indices(dofs_per_cell);
1387 dof_cell->get_mg_dof_indices(dof_indices);
1388 std::vector<double> values(dofs_per_cell);
1389 VectorHelper<VectorType>::extract(*vector,
1390 dof_indices,
1391 extract_component,
1392 values);
1393 std::fill(patch_values.begin(), patch_values.end(), 0.0);
1394
1395 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1396 for (unsigned int q_point = 0; q_point < patch_values.size(); ++q_point)
1397 patch_values[q_point] +=
1398 values[i] * fe_patch_values.shape_value(i, q_point);
1399 }
1400
1401
1402
1403 template <typename DoFHandlerType, typename VectorType>
1404 void
get_function_values(const FEValuesBase<DoFHandlerType::dimension,DoFHandlerType::space_dimension> & fe_patch_values,const ComponentExtractor extract_component,std::vector<dealii::Vector<double>> & patch_values_system)1405 MGDataEntry<DoFHandlerType, VectorType>::get_function_values(
1406 const FEValuesBase<DoFHandlerType::dimension,
1407 DoFHandlerType::space_dimension> &fe_patch_values,
1408 const ComponentExtractor extract_component,
1409 std::vector<dealii::Vector<double>> &patch_values_system) const
1410 {
1411 Assert(extract_component == ComponentExtractor::real_part,
1412 ExcNotImplemented());
1413
1414 typename DoFHandlerType::level_cell_iterator dof_cell(
1415 &fe_patch_values.get_cell()->get_triangulation(),
1416 fe_patch_values.get_cell()->level(),
1417 fe_patch_values.get_cell()->index(),
1418 this->dof_handler);
1419
1420 const VectorType *vector = &((*vectors)[dof_cell->level()]);
1421
1422 const unsigned int dofs_per_cell =
1423 this->dof_handler->get_fe(0).n_dofs_per_cell();
1424
1425 std::vector<types::global_dof_index> dof_indices(dofs_per_cell);
1426 dof_cell->get_mg_dof_indices(dof_indices);
1427 std::vector<double> values(dofs_per_cell);
1428 VectorHelper<VectorType>::extract(*vector,
1429 dof_indices,
1430 extract_component,
1431 values);
1432
1433 const unsigned int n_components = fe_patch_values.get_fe().n_components();
1434 const unsigned int n_eval_points = fe_patch_values.n_quadrature_points;
1435
1436 AssertDimension(patch_values_system.size(), n_eval_points);
1437 for (unsigned int q = 0; q < n_eval_points; q++)
1438 {
1439 AssertDimension(patch_values_system[q].size(), n_components);
1440 patch_values_system[q] = 0.0;
1441
1442 for (unsigned int i = 0; i < dofs_per_cell; ++i)
1443 for (unsigned int c = 0; c < n_components; ++c)
1444 patch_values_system[q](c) +=
1445 values[i] * fe_patch_values.shape_value_component(i, q, c);
1446 }
1447 }
1448
1449 } // namespace DataOutImplementation
1450 } // namespace internal
1451
1452
1453
1454 template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
DataOut_DoFData()1455 DataOut_DoFData<DoFHandlerType, patch_dim, patch_space_dim>::DataOut_DoFData()
1456 : triangulation(nullptr, typeid(*this).name())
1457 , dofs(nullptr, typeid(*this).name())
1458 {}
1459
1460
1461
1462 template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
~DataOut_DoFData()1463 DataOut_DoFData<DoFHandlerType, patch_dim, patch_space_dim>::~DataOut_DoFData()
1464 {
1465 // virtual functions called in constructors and destructors never use the
1466 // override in a derived class for clarity be explicit on which function is
1467 // called
1468 DataOut_DoFData<DoFHandlerType, patch_dim, patch_space_dim>::clear();
1469 }
1470
1471
1472
1473 template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
1474 void
attach_dof_handler(const DoFHandlerType & d)1475 DataOut_DoFData<DoFHandlerType, patch_dim, patch_space_dim>::attach_dof_handler(
1476 const DoFHandlerType &d)
1477 {
1478 Assert(dof_data.size() == 0,
1479 Exceptions::DataOutImplementation::ExcOldDataStillPresent());
1480 Assert(cell_data.size() == 0,
1481 Exceptions::DataOutImplementation::ExcOldDataStillPresent());
1482
1483 triangulation =
1484 SmartPointer<const Triangulation<DoFHandlerType::dimension,
1485 DoFHandlerType::space_dimension>>(
1486 &d.get_triangulation(), typeid(*this).name());
1487 dofs = SmartPointer<const DoFHandlerType>(&d, typeid(*this).name());
1488 }
1489
1490
1491
1492 template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
1493 void
1494 DataOut_DoFData<DoFHandlerType, patch_dim, patch_space_dim>::
attach_triangulation(const Triangulation<DoFHandlerType::dimension,DoFHandlerType::space_dimension> & tria)1495 attach_triangulation(
1496 const Triangulation<DoFHandlerType::dimension,
1497 DoFHandlerType::space_dimension> &tria)
1498 {
1499 Assert(dof_data.size() == 0,
1500 Exceptions::DataOutImplementation::ExcOldDataStillPresent());
1501 Assert(cell_data.size() == 0,
1502 Exceptions::DataOutImplementation::ExcOldDataStillPresent());
1503
1504 triangulation =
1505 SmartPointer<const Triangulation<DoFHandlerType::dimension,
1506 DoFHandlerType::space_dimension>>(
1507 &tria, typeid(*this).name());
1508 }
1509
1510
1511
1512 template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
1513 template <typename VectorType>
1514 void
add_data_vector(const DoFHandlerType & dof_handler,const VectorType & vec,const DataPostprocessor<DoFHandlerType::space_dimension> & data_postprocessor)1515 DataOut_DoFData<DoFHandlerType, patch_dim, patch_space_dim>::add_data_vector(
1516 const DoFHandlerType & dof_handler,
1517 const VectorType & vec,
1518 const DataPostprocessor<DoFHandlerType::space_dimension> &data_postprocessor)
1519 {
1520 // this is a specialized version of the other function where we have a
1521 // postprocessor. if we do, we know that we have type_dof_data, which makes
1522 // things a bit simpler, we also don't need to deal with some of the other
1523 // stuff and use a different constructor of DataEntry
1524 if (triangulation != nullptr)
1525 {
1526 Assert(&dof_handler.get_triangulation() == triangulation,
1527 ExcMessage("The triangulation attached to the DoFHandler does not "
1528 "match with the one set previously"));
1529 }
1530 else
1531 {
1532 triangulation =
1533 SmartPointer<const Triangulation<DoFHandlerType::dimension,
1534 DoFHandlerType::space_dimension>>(
1535 &dof_handler.get_triangulation(), typeid(*this).name());
1536 }
1537
1538 Assert(vec.size() == dof_handler.n_dofs(),
1539 Exceptions::DataOutImplementation::ExcInvalidVectorSize(
1540 vec.size(),
1541 dof_handler.n_dofs(),
1542 dof_handler.get_triangulation().n_active_cells()));
1543
1544
1545 auto new_entry = std::make_unique<
1546 internal::DataOutImplementation::DataEntry<DoFHandlerType, VectorType>>(
1547 &dof_handler, &vec, &data_postprocessor);
1548 dof_data.emplace_back(std::move(new_entry));
1549 }
1550
1551
1552
1553 template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
1554 template <typename VectorType>
1555 void
1556 DataOut_DoFData<DoFHandlerType, patch_dim, patch_space_dim>::
add_data_vector_internal(const DoFHandlerType * dof_handler,const VectorType & data_vector,const std::vector<std::string> & names,const DataVectorType type,const std::vector<DataComponentInterpretation::DataComponentInterpretation> & data_component_interpretation_,const bool deduce_output_names)1557 add_data_vector_internal(
1558 const DoFHandlerType * dof_handler,
1559 const VectorType & data_vector,
1560 const std::vector<std::string> &names,
1561 const DataVectorType type,
1562 const std::vector<DataComponentInterpretation::DataComponentInterpretation>
1563 & data_component_interpretation_,
1564 const bool deduce_output_names)
1565 {
1566 // Check available mesh information:
1567 if (triangulation == nullptr)
1568 {
1569 Assert(dof_handler != nullptr, ExcInternalError());
1570 triangulation =
1571 SmartPointer<const Triangulation<DoFHandlerType::dimension,
1572 DoFHandlerType::space_dimension>>(
1573 &dof_handler->get_triangulation(), typeid(*this).name());
1574 }
1575
1576 if (dof_handler != nullptr)
1577 {
1578 Assert(&dof_handler->get_triangulation() == triangulation,
1579 ExcMessage("The triangulation attached to the DoFHandler does not "
1580 "match with the one set previously"));
1581 }
1582
1583 // Figure out the data type:
1584 DataVectorType actual_type = type;
1585 if (type == type_automatic)
1586 {
1587 Assert(
1588 (dof_handler == nullptr) ||
1589 (triangulation->n_active_cells() != dof_handler->n_dofs()),
1590 ExcMessage(
1591 "Unable to determine the type of vector automatically because the number of DoFs "
1592 "is equal to the number of cells. Please specify DataVectorType."));
1593
1594 if (data_vector.size() == triangulation->n_active_cells())
1595 actual_type = type_cell_data;
1596 else
1597 actual_type = type_dof_data;
1598 }
1599 Assert(actual_type == type_cell_data || actual_type == type_dof_data,
1600 ExcInternalError());
1601
1602 // If necessary, append '_1', '_2', etc. to component names:
1603 std::vector<std::string> deduced_names;
1604 if (deduce_output_names && actual_type == type_dof_data)
1605 {
1606 Assert(names.size() == 1, ExcInternalError());
1607 Assert(dof_handler != nullptr, ExcInternalError());
1608 Assert(dof_handler->n_dofs() > 0,
1609 ExcMessage("The DoF handler attached to the current output vector "
1610 "does not have any degrees of freedom, so it is not "
1611 "possible to output DoF data in this context."));
1612 const std::string name = names[0];
1613 const unsigned int n_components = dof_handler->get_fe(0).n_components();
1614 deduced_names.resize(n_components);
1615 if (n_components > 1)
1616 {
1617 for (unsigned int i = 0; i < n_components; ++i)
1618 {
1619 deduced_names[i] = name + '_' + std::to_string(i);
1620 }
1621 }
1622 else
1623 {
1624 deduced_names[0] = name;
1625 }
1626 }
1627 else
1628 {
1629 deduced_names = names;
1630 }
1631
1632 // Check that things have the right sizes for the data type:
1633 switch (actual_type)
1634 {
1635 case type_cell_data:
1636 Assert(data_vector.size() == triangulation->n_active_cells(),
1637 ExcDimensionMismatch(data_vector.size(),
1638 triangulation->n_active_cells()));
1639 Assert(deduced_names.size() == 1,
1640 Exceptions::DataOutImplementation::ExcInvalidNumberOfNames(
1641 deduced_names.size(), 1));
1642 break;
1643 case type_dof_data:
1644 Assert(dof_handler != nullptr,
1645 Exceptions::DataOutImplementation::ExcNoDoFHandlerSelected());
1646 Assert(data_vector.size() == dof_handler->n_dofs(),
1647 Exceptions::DataOutImplementation::ExcInvalidVectorSize(
1648 data_vector.size(),
1649 dof_handler->n_dofs(),
1650 triangulation->n_active_cells()));
1651 Assert(deduced_names.size() == dof_handler->get_fe(0).n_components(),
1652 Exceptions::DataOutImplementation::ExcInvalidNumberOfNames(
1653 deduced_names.size(), dof_handler->get_fe(0).n_components()));
1654 break;
1655 default:
1656 Assert(false, ExcInternalError());
1657 }
1658
1659 const auto &data_component_interpretation =
1660 (data_component_interpretation_.size() != 0 ?
1661 data_component_interpretation_ :
1662 std::vector<DataComponentInterpretation::DataComponentInterpretation>(
1663 deduced_names.size(),
1664 DataComponentInterpretation::component_is_scalar));
1665
1666 // finally, add the data vector:
1667 auto new_entry = std::make_unique<
1668 internal::DataOutImplementation::DataEntry<DoFHandlerType, VectorType>>(
1669 dof_handler, &data_vector, deduced_names, data_component_interpretation);
1670
1671 if (actual_type == type_dof_data)
1672 dof_data.emplace_back(std::move(new_entry));
1673 else
1674 cell_data.emplace_back(std::move(new_entry));
1675 }
1676
1677
1678
1679 template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
1680 template <class VectorType>
1681 void
add_mg_data_vector(const DoFHandlerType & dof_handler,const MGLevelObject<VectorType> & data,const std::string & name)1682 DataOut_DoFData<DoFHandlerType, patch_dim, patch_space_dim>::add_mg_data_vector(
1683 const DoFHandlerType & dof_handler,
1684 const MGLevelObject<VectorType> &data,
1685 const std::string & name)
1686 {
1687 // forward the call to the vector version:
1688 std::vector<std::string> names(1, name);
1689 add_mg_data_vector(dof_handler, data, names);
1690 }
1691
1692
1693
1694 template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
1695 template <class VectorType>
1696 void
add_mg_data_vector(const DoFHandlerType & dof_handler,const MGLevelObject<VectorType> & data,const std::vector<std::string> & names,const std::vector<DataComponentInterpretation::DataComponentInterpretation> & data_component_interpretation_)1697 DataOut_DoFData<DoFHandlerType, patch_dim, patch_space_dim>::add_mg_data_vector(
1698 const DoFHandlerType & dof_handler,
1699 const MGLevelObject<VectorType> &data,
1700 const std::vector<std::string> & names,
1701 const std::vector<DataComponentInterpretation::DataComponentInterpretation>
1702 &data_component_interpretation_)
1703 {
1704 if (triangulation == nullptr)
1705 triangulation =
1706 SmartPointer<const Triangulation<DoFHandlerType::dimension,
1707 DoFHandlerType::space_dimension>>(
1708 &dof_handler.get_triangulation(), typeid(*this).name());
1709
1710 Assert(&dof_handler.get_triangulation() == triangulation,
1711 ExcMessage("The triangulation attached to the DoFHandler does not "
1712 "match with the one set previously"));
1713
1714 const unsigned int n_components = dof_handler.get_fe(0).n_components();
1715 std::vector<std::string> deduced_names = names;
1716
1717 if (names.size() == 1 && n_components > 1)
1718 {
1719 deduced_names.resize(n_components);
1720 for (unsigned int i = 0; i < n_components; ++i)
1721 {
1722 deduced_names[i] = names[0] + '_' + std::to_string(i);
1723 }
1724 }
1725
1726 Assert(deduced_names.size() == n_components,
1727 ExcMessage("Invalid number of names given."));
1728
1729 const std::vector<DataComponentInterpretation::DataComponentInterpretation>
1730 &data_component_interpretation =
1731 (data_component_interpretation_.size() != 0 ?
1732 data_component_interpretation_ :
1733 std::vector<DataComponentInterpretation::DataComponentInterpretation>(
1734 n_components, DataComponentInterpretation::component_is_scalar));
1735
1736 Assert(data_component_interpretation.size() == n_components,
1737 ExcMessage(
1738 "Invalid number of entries in data_component_interpretation."));
1739
1740 auto new_entry = std::make_unique<
1741 internal::DataOutImplementation::MGDataEntry<DoFHandlerType, VectorType>>(
1742 &dof_handler, &data, deduced_names, data_component_interpretation);
1743 dof_data.emplace_back(std::move(new_entry));
1744 }
1745
1746
1747
1748 template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
1749 void
1750 DataOut_DoFData<DoFHandlerType, patch_dim, patch_space_dim>::
clear_data_vectors()1751 clear_data_vectors()
1752 {
1753 dof_data.erase(dof_data.begin(), dof_data.end());
1754 cell_data.erase(cell_data.begin(), cell_data.end());
1755
1756 // delete patches
1757 std::vector<Patch> dummy;
1758 patches.swap(dummy);
1759 }
1760
1761
1762
1763 template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
1764 void
1765 DataOut_DoFData<DoFHandlerType, patch_dim, patch_space_dim>::
clear_input_data_references()1766 clear_input_data_references()
1767 {
1768 for (unsigned int i = 0; i < dof_data.size(); ++i)
1769 dof_data[i]->clear();
1770
1771 for (unsigned int i = 0; i < cell_data.size(); ++i)
1772 cell_data[i]->clear();
1773
1774 if (dofs != nullptr)
1775 dofs = nullptr;
1776 }
1777
1778
1779
1780 template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
1781 void
clear()1782 DataOut_DoFData<DoFHandlerType, patch_dim, patch_space_dim>::clear()
1783 {
1784 dof_data.erase(dof_data.begin(), dof_data.end());
1785 cell_data.erase(cell_data.begin(), cell_data.end());
1786
1787 if (dofs != nullptr)
1788 dofs = nullptr;
1789
1790 // delete patches
1791 std::vector<Patch> dummy;
1792 patches.swap(dummy);
1793 }
1794
1795
1796
1797 template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
1798 std::vector<std::string>
get_dataset_names()1799 DataOut_DoFData<DoFHandlerType, patch_dim, patch_space_dim>::get_dataset_names()
1800 const
1801 {
1802 std::vector<std::string> names;
1803
1804 // Loop over all DoF-data datasets and push the names. If the
1805 // vector underlying a data set is complex-valued, then
1806 // expand it into its real and imaginary part. Note, however,
1807 // that what comes back from a postprocessor is *always*
1808 // real-valued, regardless of what goes in, so we don't
1809 // have this to do this name expansion for data sets that
1810 // have a postprocessor.
1811 //
1812 // The situation is made complicated when considering vector- and
1813 // tensor-valued component sets. This is because if, for example, we have a
1814 // complex-valued vector, we don't want to output Re(u_x), then Im(u_x), then
1815 // Re(u_y), etc. That's because if we did this, then visualization programs
1816 // will not easily be able to patch together the 1st, 3rd, 5th components into
1817 // the vector representing the real part of a vector field, and similarly for
1818 // the 2nd, 4th, 6th component for the imaginary part of the vector field.
1819 // Rather, we need to put all real components of the same vector field into
1820 // consecutive components.
1821 //
1822 // This sort of logic is also explained in some detail in
1823 // DataOut::build_one_patch().
1824 for (const auto &input_data : dof_data)
1825 if (input_data->is_complex_valued() == false ||
1826 (input_data->postprocessor != nullptr))
1827 {
1828 for (const auto &name : input_data->names)
1829 names.push_back(name);
1830 }
1831 else
1832 {
1833 // OK, so we have a complex-valued vector. We then need to go through
1834 // all components and order them appropriately
1835 for (unsigned int i = 0; i < input_data->names.size();
1836 /* increment of i happens below */)
1837 {
1838 switch (input_data->data_component_interpretation[i])
1839 {
1840 case DataComponentInterpretation::component_is_scalar:
1841 {
1842 // It's a scalar. Just output real and imaginary parts one
1843 // after the other:
1844 names.push_back(input_data->names[i] + "_re");
1845 names.push_back(input_data->names[i] + "_im");
1846
1847 // Move forward by one component
1848 ++i;
1849
1850 break;
1851 }
1852
1853 case DataComponentInterpretation::component_is_part_of_vector:
1854 {
1855 // It's a vector. First output all real parts, then all
1856 // imaginary parts:
1857 const unsigned int size = patch_space_dim;
1858 for (unsigned int vec_comp = 0; vec_comp < size; ++vec_comp)
1859 names.push_back(input_data->names[i + vec_comp] + "_re");
1860 for (unsigned int vec_comp = 0; vec_comp < size; ++vec_comp)
1861 names.push_back(input_data->names[i + vec_comp] + "_im");
1862
1863 // Move forward by dim components
1864 i += size;
1865
1866 break;
1867 }
1868
1869 case DataComponentInterpretation::component_is_part_of_tensor:
1870 {
1871 // It's a tensor. First output all real parts, then all
1872 // imaginary parts:
1873 const unsigned int size = patch_space_dim * patch_space_dim;
1874 for (unsigned int tensor_comp = 0; tensor_comp < size;
1875 ++tensor_comp)
1876 names.push_back(input_data->names[i + tensor_comp] +
1877 "_re");
1878 for (unsigned int tensor_comp = 0; tensor_comp < size;
1879 ++tensor_comp)
1880 names.push_back(input_data->names[i + tensor_comp] +
1881 "_im");
1882
1883 // Move forward by dim*dim components
1884 i += size;
1885
1886 break;
1887 }
1888
1889 default:
1890 Assert(false, ExcInternalError());
1891 }
1892 }
1893 }
1894
1895 // Do the same as above for cell-type data. This is simpler because it
1896 // is always scalar, and so we don't have to worry about whether some
1897 // components together form vectors or tensors.
1898 for (const auto &input_data : cell_data)
1899 {
1900 Assert(input_data->names.size() == 1, ExcInternalError());
1901 if ((input_data->is_complex_valued() == false) ||
1902 (input_data->postprocessor != nullptr))
1903 names.push_back(input_data->names[0]);
1904 else
1905 {
1906 names.push_back(input_data->names[0] + "_re");
1907 names.push_back(input_data->names[0] + "_im");
1908 }
1909 }
1910
1911 return names;
1912 }
1913
1914
1915
1916 template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
1917 std::vector<
1918 std::tuple<unsigned int,
1919 unsigned int,
1920 std::string,
1921 DataComponentInterpretation::DataComponentInterpretation>>
1922 DataOut_DoFData<DoFHandlerType, patch_dim, patch_space_dim>::
get_nonscalar_data_ranges()1923 get_nonscalar_data_ranges() const
1924 {
1925 std::vector<
1926 std::tuple<unsigned int,
1927 unsigned int,
1928 std::string,
1929 DataComponentInterpretation::DataComponentInterpretation>>
1930 ranges;
1931
1932 // collect the ranges of dof and cell data
1933 unsigned int output_component = 0;
1934 for (const auto &input_data : dof_data)
1935 for (unsigned int i = 0; i < input_data->n_output_variables;
1936 /* i is updated below */)
1937 // see what kind of data we have here. note that for the purpose of the
1938 // current function all we care about is vector data
1939 switch (input_data->data_component_interpretation[i])
1940 {
1941 case DataComponentInterpretation::component_is_scalar:
1942 {
1943 // Just move component forward by one (or two if the
1944 // component happens to be complex-valued and we don't use a
1945 // postprocessor)
1946 // -- postprocessors always return real-valued things)
1947 ++i;
1948 output_component += (input_data->is_complex_valued() &&
1949 (input_data->postprocessor == nullptr) ?
1950 2 :
1951 1);
1952
1953 break;
1954 }
1955
1956 case DataComponentInterpretation::component_is_part_of_vector:
1957 {
1958 // ensure that there is a continuous number of next space_dim
1959 // components that all deal with vectors
1960 Assert(
1961 i + patch_space_dim <= input_data->n_output_variables,
1962 Exceptions::DataOutImplementation::ExcInvalidVectorDeclaration(
1963 i, input_data->names[i]));
1964 for (unsigned int dd = 1; dd < patch_space_dim; ++dd)
1965 Assert(
1966 input_data->data_component_interpretation[i + dd] ==
1967 DataComponentInterpretation::component_is_part_of_vector,
1968 Exceptions::DataOutImplementation::
1969 ExcInvalidVectorDeclaration(i, input_data->names[i]));
1970
1971 // all seems right, so figure out whether there is a common
1972 // name to these components. if not, leave the name empty and
1973 // let the output format writer decide what to do here
1974 std::string name = input_data->names[i];
1975 for (unsigned int dd = 1; dd < patch_space_dim; ++dd)
1976 if (name != input_data->names[i + dd])
1977 {
1978 name = "";
1979 break;
1980 }
1981
1982 // Finally add a corresponding range. If this is a real-valued
1983 // vector, then we only need to do this once. But if it is a
1984 // complex-valued vector and it is not postprocessed, then we need
1985 // to do it twice -- once for the real parts and once for the
1986 // imaginary parts
1987 //
1988 // This sort of logic is also explained in some detail in
1989 // DataOut::build_one_patch().
1990 if (input_data->is_complex_valued() == false ||
1991 (input_data->postprocessor != nullptr))
1992 {
1993 ranges.emplace_back(std::forward_as_tuple(
1994 output_component,
1995 output_component + patch_space_dim - 1,
1996 name,
1997 DataComponentInterpretation::component_is_part_of_vector));
1998
1999 // increase the 'component' counter by the appropriate amount,
2000 // same for 'i', since we have already dealt with all these
2001 // components
2002 output_component += patch_space_dim;
2003 i += patch_space_dim;
2004 }
2005 else
2006 {
2007 ranges.emplace_back(std::forward_as_tuple(
2008 output_component,
2009 output_component + patch_space_dim - 1,
2010 name + "_re",
2011 DataComponentInterpretation::component_is_part_of_vector));
2012 output_component += patch_space_dim;
2013
2014 ranges.emplace_back(std::forward_as_tuple(
2015 output_component,
2016 output_component + patch_space_dim - 1,
2017 name + "_im",
2018 DataComponentInterpretation::component_is_part_of_vector));
2019 output_component += patch_space_dim;
2020
2021 i += patch_space_dim;
2022 }
2023
2024
2025 break;
2026 }
2027
2028 case DataComponentInterpretation::component_is_part_of_tensor:
2029 {
2030 const unsigned int size = patch_space_dim * patch_space_dim;
2031 // ensure that there is a continuous number of next
2032 // space_dim*space_dim components that all deal with tensors
2033 Assert(
2034 i + size <= input_data->n_output_variables,
2035 Exceptions::DataOutImplementation::ExcInvalidTensorDeclaration(
2036 i, input_data->names[i]));
2037 for (unsigned int dd = 1; dd < size; ++dd)
2038 Assert(
2039 input_data->data_component_interpretation[i + dd] ==
2040 DataComponentInterpretation::component_is_part_of_tensor,
2041 Exceptions::DataOutImplementation::
2042 ExcInvalidTensorDeclaration(i, input_data->names[i]));
2043
2044 // all seems right, so figure out whether there is a common
2045 // name to these components. if not, leave the name empty and
2046 // let the output format writer decide what to do here
2047 std::string name = input_data->names[i];
2048 for (unsigned int dd = 1; dd < size; ++dd)
2049 if (name != input_data->names[i + dd])
2050 {
2051 name = "";
2052 break;
2053 }
2054
2055 // Finally add a corresponding range. If this is a real-valued
2056 // tensor, then we only need to do this once. But if it is a
2057 // complex-valued tensor and it is not postprocessed, then we need
2058 // to do it twice -- once for the real parts and once for the
2059 // imaginary parts
2060 //
2061 // This sort of logic is also explained in some detail in
2062 // DataOut::build_one_patch().
2063 if (input_data->is_complex_valued() == false ||
2064 (input_data->postprocessor != nullptr))
2065 {
2066 ranges.emplace_back(std::forward_as_tuple(
2067 output_component,
2068 output_component + size - 1,
2069 name,
2070 DataComponentInterpretation::component_is_part_of_tensor));
2071
2072 // increase the 'component' counter by the appropriate amount,
2073 // same for 'i', since we have already dealt with all these
2074 // components
2075 output_component += size;
2076 i += size;
2077 }
2078 else
2079 {
2080 ranges.emplace_back(std::forward_as_tuple(
2081 output_component,
2082 output_component + size - 1,
2083 name + "_re",
2084 DataComponentInterpretation::component_is_part_of_tensor));
2085 output_component += size;
2086
2087 ranges.emplace_back(std::forward_as_tuple(
2088 output_component,
2089 output_component + size - 1,
2090 name + "_im",
2091 DataComponentInterpretation::component_is_part_of_tensor));
2092 output_component += size;
2093
2094 i += size;
2095 }
2096 break;
2097 }
2098
2099 default:
2100 Assert(false, ExcNotImplemented());
2101 }
2102
2103 // note that we do not have to traverse the list of cell data here because
2104 // cell data is one value per (logical) cell and therefore cannot be a
2105 // vector
2106
2107 return ranges;
2108 }
2109
2110
2111
2112 template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
2113 const std::vector<dealii::DataOutBase::Patch<patch_dim, patch_space_dim>> &
get_patches()2114 DataOut_DoFData<DoFHandlerType, patch_dim, patch_space_dim>::get_patches() const
2115 {
2116 return patches;
2117 }
2118
2119
2120
2121 template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
2122 std::vector<
2123 std::shared_ptr<dealii::hp::FECollection<DoFHandlerType::dimension,
2124 DoFHandlerType::space_dimension>>>
get_fes()2125 DataOut_DoFData<DoFHandlerType, patch_dim, patch_space_dim>::get_fes() const
2126 {
2127 const unsigned int dhdim = DoFHandlerType::dimension;
2128 const unsigned int dhspacedim = DoFHandlerType::space_dimension;
2129 std::vector<std::shared_ptr<dealii::hp::FECollection<dhdim, dhspacedim>>>
2130 finite_elements(this->dof_data.size());
2131 for (unsigned int i = 0; i < this->dof_data.size(); ++i)
2132 {
2133 Assert(dof_data[i]->dof_handler != nullptr,
2134 Exceptions::DataOutImplementation::ExcNoDoFHandlerSelected());
2135
2136 // avoid creating too many finite elements and doing a lot of work on
2137 // initializing FEValues downstream: if two DoFHandlers are the same
2138 // (checked by pointer comparison), we can re-use the shared_ptr object
2139 // for the second one. We cannot check for finite element equalities
2140 // because we need different FEValues objects for different dof
2141 // handlers.
2142 bool duplicate = false;
2143 for (unsigned int j = 0; j < i; ++j)
2144 if (dof_data[i]->dof_handler == dof_data[j]->dof_handler)
2145 {
2146 finite_elements[i] = finite_elements[j];
2147 duplicate = true;
2148 }
2149 if (duplicate == false)
2150 finite_elements[i] =
2151 std::make_shared<dealii::hp::FECollection<dhdim, dhspacedim>>(
2152 this->dof_data[i]->dof_handler->get_fe_collection());
2153 }
2154 if (this->dof_data.empty())
2155 {
2156 finite_elements.resize(1);
2157 finite_elements[0] =
2158 std::make_shared<dealii::hp::FECollection<dhdim, dhspacedim>>(
2159 FE_DGQ<dhdim, dhspacedim>(0));
2160 }
2161 return finite_elements;
2162 }
2163
2164
2165
2166 template <typename DoFHandlerType, int patch_dim, int patch_space_dim>
2167 std::size_t
2168 DataOut_DoFData<DoFHandlerType, patch_dim, patch_space_dim>::
memory_consumption()2169 memory_consumption() const
2170 {
2171 return (DataOutInterface<patch_dim, patch_space_dim>::memory_consumption() +
2172 MemoryConsumption::memory_consumption(dofs) +
2173 MemoryConsumption::memory_consumption(patches));
2174 }
2175
2176 DEAL_II_NAMESPACE_CLOSE
2177
2178 #endif
2179