Reference documentation for deal.II version 9.3.3
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
data_out.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1999 - 2021 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
17
20
21#include <deal.II/fe/fe.h>
22#include <deal.II/fe/fe_dgq.h>
25
26#include <deal.II/grid/tria.h>
28
31
33
34#include <sstream>
35
37
38
39namespace internal
40{
41 namespace DataOutImplementation
42 {
43 template <int dim, int spacedim>
45 const unsigned int n_datasets,
46 const unsigned int n_subdivisions,
47 const std::vector<unsigned int> &n_postprocessor_outputs,
48 const ::hp::MappingCollection<dim, spacedim> &mapping,
49 const std::vector<
50 std::shared_ptr<::hp::FECollection<dim, spacedim>>>
51 & finite_elements,
52 const UpdateFlags update_flags,
53 const std::vector<std::vector<unsigned int>> &cell_to_patch_index_map)
54 : ParallelDataBase<dim, spacedim>(n_datasets,
55 n_subdivisions,
56 n_postprocessor_outputs,
57 mapping,
58 finite_elements,
59 update_flags,
60 false)
61 , cell_to_patch_index_map(&cell_to_patch_index_map)
62 {}
63 } // namespace DataOutImplementation
64} // namespace internal
65
66
67
68template <int dim, typename DoFHandlerType>
70{
71 // For the moment, just call the existing virtual functions. This
72 // preserves backward compatibility. When these deprecated functions are
73 // removed, we'll just inline their functionality into the lambda
74 // functions created here.
75 set_cell_selection(
76 [this](const Triangulation<dim, spacedim> &) {
77 return this->first_locally_owned_cell();
78 },
79 [this](const Triangulation<dim, spacedim> &, const cell_iterator &cell) {
80 return this->next_locally_owned_cell(cell);
81 });
82}
83
84
85
86template <int dim, typename DoFHandlerType>
87void
89 const std::pair<cell_iterator, unsigned int> *cell_and_index,
90 internal::DataOutImplementation::ParallelData<DoFHandlerType::dimension,
91 DoFHandlerType::space_dimension>
92 & scratch_data,
93 const unsigned int n_subdivisions,
94 const CurvedCellRegion curved_cell_region)
95{
96 // first create the output object that we will write into
97 ::DataOutBase::Patch<DoFHandlerType::dimension,
98 DoFHandlerType::space_dimension>
99 patch;
100 patch.n_subdivisions = n_subdivisions;
101 patch.reference_cell = cell_and_index->first->reference_cell();
102
103 // initialize FEValues
104 scratch_data.reinit_all_fe_values(this->dof_data, cell_and_index->first);
105
107 &fe_patch_values = scratch_data.get_present_fe_values(0);
108
109 // set the vertices of the patch. if the mapping does not preserve locations
110 // (e.g. MappingQEulerian), we need to compute the offset of the vertex for
111 // the graphical output. Otherwise, we can just use the vertex info.
112 for (const unsigned int vertex : cell_and_index->first->vertex_indices())
113 if (scratch_data.mapping_collection[0].preserves_vertex_locations())
114 patch.vertices[vertex] = cell_and_index->first->vertex(vertex);
115 else
116 patch.vertices[vertex] =
118 cell_and_index->first,
120
121 const unsigned int n_q_points = fe_patch_values.n_quadrature_points;
122
123 scratch_data.patch_values_scalar.solution_values.resize(n_q_points);
124 scratch_data.patch_values_scalar.solution_gradients.resize(n_q_points);
125 scratch_data.patch_values_scalar.solution_hessians.resize(n_q_points);
126 scratch_data.patch_values_system.solution_values.resize(n_q_points);
127 scratch_data.patch_values_system.solution_gradients.resize(n_q_points);
128 scratch_data.patch_values_system.solution_hessians.resize(n_q_points);
129
130 for (unsigned int dataset = 0;
131 dataset < scratch_data.postprocessed_values.size();
132 ++dataset)
133 if (scratch_data.postprocessed_values[dataset].size() != 0)
134 scratch_data.postprocessed_values[dataset].resize(n_q_points);
135
136 // First fill the geometric information for the patch: Where are the
137 // nodes in question located.
138 //
139 // Depending on the requested output of curved cells, if necessary
140 // append the quadrature points to the last rows of the patch.data
141 // member. This is the case if we want to produce curved cells at the
142 // boundary and this cell actually is at the boundary, or else if we
143 // want to produce curved cells everywhere
144 //
145 // note: a cell is *always* at the boundary if dim<spacedim
146 if (curved_cell_region == curved_inner_cells ||
147 (curved_cell_region == curved_boundary &&
148 (cell_and_index->first->at_boundary() ||
149 (DoFHandlerType::dimension != DoFHandlerType::space_dimension))) ||
150 (cell_and_index->first->reference_cell() !=
151 ReferenceCells::get_hypercube<dim>()))
152 {
153 Assert(patch.space_dim == DoFHandlerType::space_dimension,
155
156 // set the flag indicating that for this cell the points are
157 // explicitly given
158 patch.points_are_available = true;
159
160 // then resize the patch.data member in order to have enough memory for
161 // the quadrature points as well, and copy the quadrature points there
162 const std::vector<Point<DoFHandlerType::space_dimension>> &q_points =
163 fe_patch_values.get_quadrature_points();
164
165 patch.data.reinit(scratch_data.n_datasets +
166 DoFHandlerType::space_dimension,
167 n_q_points);
168 for (unsigned int i = 0; i < DoFHandlerType::space_dimension; ++i)
169 for (unsigned int q = 0; q < n_q_points; ++q)
170 patch.data(patch.data.size(0) - DoFHandlerType::space_dimension + i,
171 q) = q_points[q][i];
172 }
173 else
174 {
175 patch.data.reinit(scratch_data.n_datasets, n_q_points);
176 patch.points_are_available = false;
177 }
178
179
180 // Next fill the information we get from DoF data
181 if (scratch_data.n_datasets > 0)
182 {
183 // counter for data records
184 unsigned int offset = 0;
185
186 // first fill dof_data
187 unsigned int dataset_number = 0;
188 for (const auto &dataset : this->dof_data)
189 {
190 const FEValuesBase<DoFHandlerType::dimension,
191 DoFHandlerType::space_dimension>
192 &this_fe_patch_values =
193 scratch_data.get_present_fe_values(dataset_number);
194 const unsigned int n_components =
195 this_fe_patch_values.get_fe().n_components();
196
198 *postprocessor = dataset->postprocessor;
199
200 if (postprocessor != nullptr)
201 {
202 // We have to postprocess the data, so determine, which fields
203 // have to be updated
204 const UpdateFlags update_flags =
205 postprocessor->get_needed_update_flags();
206
207 if ((n_components == 1) &&
208 (dataset->is_complex_valued() == false))
209 {
210 // At each point there is only one component of value,
211 // gradient etc. Based on the 'if' statement above, we
212 // know that the solution is scalar and real-valued, so
213 // we do not need to worry about getting any imaginary
214 // components to the postprocessor, and we can safely
215 // call the function that evaluates a scalar field
216 if (update_flags & update_values)
217 dataset->get_function_values(
218 this_fe_patch_values,
220 real_part,
221 scratch_data.patch_values_scalar.solution_values);
222
223 if (update_flags & update_gradients)
224 dataset->get_function_gradients(
225 this_fe_patch_values,
227 real_part,
228 scratch_data.patch_values_scalar.solution_gradients);
229
230 if (update_flags & update_hessians)
231 dataset->get_function_hessians(
232 this_fe_patch_values,
234 real_part,
235 scratch_data.patch_values_scalar.solution_hessians);
236
237 // Also fill some of the other fields postprocessors may
238 // want to access.
239 if (update_flags & update_quadrature_points)
240 scratch_data.patch_values_scalar.evaluation_points =
241 this_fe_patch_values.get_quadrature_points();
242
243 const typename DoFHandlerType::cell_iterator dh_cell(
244 &cell_and_index->first->get_triangulation(),
245 cell_and_index->first->level(),
246 cell_and_index->first->index(),
247 dataset->dof_handler);
248 scratch_data.patch_values_scalar.template set_cell<dim>(
249 dh_cell);
250
251 // Finally call the postprocessor's function that
252 // deals with scalar inputs.
253 postprocessor->evaluate_scalar_field(
254 scratch_data.patch_values_scalar,
255 scratch_data.postprocessed_values[dataset_number]);
256 }
257 else
258 {
259 // At each point we now have to evaluate a vector valued
260 // function and its derivatives. It may be that the solution
261 // is scalar and complex-valued, but we treat this as a vector
262 // field with two components.
263
264 // At this point, we need to ask how we fill the fields that
265 // we want to pass on to the postprocessor. If the field in
266 // question is real-valued, we'll just extract the (only)
267 // real component from the solution fields
268 if (dataset->is_complex_valued() == false)
269 {
270 scratch_data.resize_system_vectors(n_components);
271
272 if (update_flags & update_values)
273 dataset->get_function_values(
274 this_fe_patch_values,
276 real_part,
277 scratch_data.patch_values_system.solution_values);
278
279 if (update_flags & update_gradients)
280 dataset->get_function_gradients(
281 this_fe_patch_values,
283 real_part,
284 scratch_data.patch_values_system.solution_gradients);
285
286 if (update_flags & update_hessians)
287 dataset->get_function_hessians(
288 this_fe_patch_values,
290 real_part,
291 scratch_data.patch_values_system.solution_hessians);
292 }
293 else
294 {
295 // The solution is complex-valued. Let's cover the scalar
296 // case first (i.e., one scalar but complex-valued field,
297 // which we will have to split into its real and imaginar
298 // parts).
299 if (n_components == 1)
300 {
301 scratch_data.resize_system_vectors(2);
302
303 // First get the real component of the scalar solution
304 // and copy the data into the
305 // scratch_data.patch_values_system output fields
306 if (update_flags & update_values)
307 {
308 dataset->get_function_values(
309 this_fe_patch_values,
311 ComponentExtractor::real_part,
312 scratch_data.patch_values_scalar
313 .solution_values);
314
315 for (unsigned int i = 0;
316 i < scratch_data.patch_values_scalar
317 .solution_values.size();
318 ++i)
319 {
321 scratch_data.patch_values_system
322 .solution_values[i]
323 .size(),
324 2);
325 scratch_data.patch_values_system
326 .solution_values[i][0] =
327 scratch_data.patch_values_scalar
328 .solution_values[i];
329 }
330 }
331
332 if (update_flags & update_gradients)
333 {
334 dataset->get_function_gradients(
335 this_fe_patch_values,
337 ComponentExtractor::real_part,
338 scratch_data.patch_values_scalar
339 .solution_gradients);
340
341 for (unsigned int i = 0;
342 i < scratch_data.patch_values_scalar
343 .solution_gradients.size();
344 ++i)
345 {
347 scratch_data.patch_values_system
348 .solution_values[i]
349 .size(),
350 2);
351 scratch_data.patch_values_system
352 .solution_gradients[i][0] =
353 scratch_data.patch_values_scalar
354 .solution_gradients[i];
355 }
356 }
357
358 if (update_flags & update_hessians)
359 {
360 dataset->get_function_hessians(
361 this_fe_patch_values,
363 ComponentExtractor::real_part,
364 scratch_data.patch_values_scalar
365 .solution_hessians);
366
367 for (unsigned int i = 0;
368 i < scratch_data.patch_values_scalar
369 .solution_hessians.size();
370 ++i)
371 {
373 scratch_data.patch_values_system
374 .solution_hessians[i]
375 .size(),
376 2);
377 scratch_data.patch_values_system
378 .solution_hessians[i][0] =
379 scratch_data.patch_values_scalar
380 .solution_hessians[i];
381 }
382 }
383
384 // Now we also have to get the imaginary
385 // component of the scalar solution
386 // and copy the data into the
387 // scratch_data.patch_values_system output fields
388 // that follow the real one
389 if (update_flags & update_values)
390 {
391 dataset->get_function_values(
392 this_fe_patch_values,
394 ComponentExtractor::imaginary_part,
395 scratch_data.patch_values_scalar
396 .solution_values);
397
398 for (unsigned int i = 0;
399 i < scratch_data.patch_values_scalar
400 .solution_values.size();
401 ++i)
402 {
404 scratch_data.patch_values_system
405 .solution_values[i]
406 .size(),
407 2);
408 scratch_data.patch_values_system
409 .solution_values[i][1] =
410 scratch_data.patch_values_scalar
411 .solution_values[i];
412 }
413 }
414
415 if (update_flags & update_gradients)
416 {
417 dataset->get_function_gradients(
418 this_fe_patch_values,
420 ComponentExtractor::imaginary_part,
421 scratch_data.patch_values_scalar
422 .solution_gradients);
423
424 for (unsigned int i = 0;
425 i < scratch_data.patch_values_scalar
426 .solution_gradients.size();
427 ++i)
428 {
430 scratch_data.patch_values_system
431 .solution_values[i]
432 .size(),
433 2);
434 scratch_data.patch_values_system
435 .solution_gradients[i][1] =
436 scratch_data.patch_values_scalar
437 .solution_gradients[i];
438 }
439 }
440
441 if (update_flags & update_hessians)
442 {
443 dataset->get_function_hessians(
444 this_fe_patch_values,
446 ComponentExtractor::imaginary_part,
447 scratch_data.patch_values_scalar
448 .solution_hessians);
449
450 for (unsigned int i = 0;
451 i < scratch_data.patch_values_scalar
452 .solution_hessians.size();
453 ++i)
454 {
456 scratch_data.patch_values_system
457 .solution_hessians[i]
458 .size(),
459 2);
460 scratch_data.patch_values_system
461 .solution_hessians[i][1] =
462 scratch_data.patch_values_scalar
463 .solution_hessians[i];
464 }
465 }
466 }
467 else
468 {
469 scratch_data.resize_system_vectors(2 * n_components);
470
471 // This is the vector-valued, complex-valued case. In
472 // essence, we just need to do the same as above,
473 // i.e., call the functions in dataset
474 // to retrieve first the real and then the imaginary
475 // part of the solution, then copy them to the
476 // scratch_data.patch_values_system. The difference to
477 // the scalar case is that there, we could (ab)use the
478 // scratch_data.patch_values_scalar members for first
479 // the real part and then the imaginary part, copying
480 // them into the scratch_data.patch_values_system
481 // variable one after the other. We can't do this
482 // here because the solution is vector-valued, and
483 // so using the single
484 // scratch_data.patch_values_system object doesn't
485 // work.
486 //
487 // Rather, we need to come up with a temporary object
488 // for this (one for each the values, the gradients,
489 // and the hessians).
490 //
491 // Compared to the previous code path, we here
492 // first get the real and imaginary parts of the
493 // values, then the real and imaginary parts of the
494 // gradients, etc. This allows us to scope the
495 // temporary objects better
496 if (update_flags & update_values)
497 {
498 std::vector<Vector<double>> tmp(
499 scratch_data.patch_values_system.solution_values
500 .size(),
501 Vector<double>(n_components));
502
503 // First get the real part into the tmp object
504 dataset->get_function_values(
505 this_fe_patch_values,
507 ComponentExtractor::real_part,
508 tmp);
509
510 // Then copy these values into the first
511 // n_components slots of the output object.
512 for (unsigned int i = 0;
513 i < scratch_data.patch_values_system
514 .solution_values.size();
515 ++i)
516 {
518 scratch_data.patch_values_system
519 .solution_values[i]
520 .size(),
521 2 * n_components);
522 for (unsigned int j = 0; j < n_components;
523 ++j)
524 scratch_data.patch_values_system
525 .solution_values[i][j] = tmp[i][j];
526 }
527
528 // Then do the same with the imaginary part,
529 // copying past the end of the previous set of
530 // values.
531 dataset->get_function_values(
532 this_fe_patch_values,
534 ComponentExtractor::imaginary_part,
535 tmp);
536
537 for (unsigned int i = 0;
538 i < scratch_data.patch_values_system
539 .solution_values.size();
540 ++i)
541 {
542 for (unsigned int j = 0; j < n_components;
543 ++j)
544 scratch_data.patch_values_system
545 .solution_values[i][j + n_components] =
546 tmp[i][j];
547 }
548 }
549
550 // Now do the exact same thing for the gradients
551 if (update_flags & update_gradients)
552 {
553 std::vector<std::vector<
555 tmp(
556 scratch_data.patch_values_system
557 .solution_gradients.size(),
558 std::vector<
560 n_components));
561
562 // First the real part
563 dataset->get_function_gradients(
564 this_fe_patch_values,
566 ComponentExtractor::real_part,
567 tmp);
568
569 for (unsigned int i = 0;
570 i < scratch_data.patch_values_system
571 .solution_gradients.size();
572 ++i)
573 {
575 scratch_data.patch_values_system
576 .solution_gradients[i]
577 .size(),
578 2 * n_components);
579 for (unsigned int j = 0; j < n_components;
580 ++j)
581 scratch_data.patch_values_system
582 .solution_gradients[i][j] = tmp[i][j];
583 }
584
585 // Then the imaginary part
586 dataset->get_function_gradients(
587 this_fe_patch_values,
589 ComponentExtractor::imaginary_part,
590 tmp);
591
592 for (unsigned int i = 0;
593 i < scratch_data.patch_values_system
594 .solution_gradients.size();
595 ++i)
596 {
597 for (unsigned int j = 0; j < n_components;
598 ++j)
599 scratch_data.patch_values_system
600 .solution_gradients[i][j + n_components] =
601 tmp[i][j];
602 }
603 }
604
605 // And finally the Hessians. Same scheme as above.
606 if (update_flags & update_hessians)
607 {
608 std::vector<std::vector<
610 tmp(
611 scratch_data.patch_values_system
612 .solution_gradients.size(),
613 std::vector<
615 n_components));
616
617 // First the real part
618 dataset->get_function_hessians(
619 this_fe_patch_values,
621 ComponentExtractor::real_part,
622 tmp);
623
624 for (unsigned int i = 0;
625 i < scratch_data.patch_values_system
626 .solution_hessians.size();
627 ++i)
628 {
630 scratch_data.patch_values_system
631 .solution_hessians[i]
632 .size(),
633 2 * n_components);
634 for (unsigned int j = 0; j < n_components;
635 ++j)
636 scratch_data.patch_values_system
637 .solution_hessians[i][j] = tmp[i][j];
638 }
639
640 // Then the imaginary part
641 dataset->get_function_hessians(
642 this_fe_patch_values,
644 ComponentExtractor::imaginary_part,
645 tmp);
646
647 for (unsigned int i = 0;
648 i < scratch_data.patch_values_system
649 .solution_hessians.size();
650 ++i)
651 {
652 for (unsigned int j = 0; j < n_components;
653 ++j)
654 scratch_data.patch_values_system
655 .solution_hessians[i][j + n_components] =
656 tmp[i][j];
657 }
658 }
659 }
660 }
661
662 // Now set other fields we may need
663 if (update_flags & update_quadrature_points)
664 scratch_data.patch_values_system.evaluation_points =
665 this_fe_patch_values.get_quadrature_points();
666
667 const typename DoFHandlerType::cell_iterator dh_cell(
668 &cell_and_index->first->get_triangulation(),
669 cell_and_index->first->level(),
670 cell_and_index->first->index(),
671 dataset->dof_handler);
672 scratch_data.patch_values_system.template set_cell<dim>(
673 dh_cell);
674
675 // Whether the solution was complex-scalar or
676 // complex-vector-valued doesn't matter -- we took it apart
677 // into several fields and so we have to call the
678 // evaluate_vector_field() function.
679 postprocessor->evaluate_vector_field(
680 scratch_data.patch_values_system,
681 scratch_data.postprocessed_values[dataset_number]);
682 }
683
684 // Now we need to copy the result of the postprocessor to
685 // the Patch object where it can then be further processed
686 // by the functions in DataOutBase
687 for (unsigned int q = 0; q < n_q_points; ++q)
688 for (unsigned int component = 0;
689 component < dataset->n_output_variables;
690 ++component)
691 patch.data(offset + component, q) =
692 scratch_data.postprocessed_values[dataset_number][q](
693 component);
694
695 // Move the counter for the output location forward as
696 // appropriate
697 offset += dataset->n_output_variables;
698 }
699 else
700 {
701 // use the given data vector directly, without a postprocessor.
702 // again, we treat single component functions separately for
703 // efficiency reasons.
704 if (n_components == 1)
705 {
706 Assert(dataset->n_output_variables == 1, ExcInternalError());
707
708 // First output the real part of the solution vector
709 dataset->get_function_values(
710 this_fe_patch_values,
712 real_part,
713 scratch_data.patch_values_scalar.solution_values);
714 for (unsigned int q = 0; q < n_q_points; ++q)
715 patch.data(offset, q) =
716 scratch_data.patch_values_scalar.solution_values[q];
717 offset += 1;
718
719 // And if there is one, also output the imaginary part. Note
720 // that the problem is scalar-valued, so we can freely add the
721 // imaginary part after the real part without having to worry
722 // that we are interleaving the real components of a vector
723 // with the imaginary components of the same vector.
724 if (dataset->is_complex_valued() == true)
725 {
726 dataset->get_function_values(
727 this_fe_patch_values,
729 imaginary_part,
730 scratch_data.patch_values_scalar.solution_values);
731 for (unsigned int q = 0; q < n_q_points; ++q)
732 patch.data(offset, q) =
733 scratch_data.patch_values_scalar.solution_values[q];
734 offset += 1;
735 }
736 }
737 else
738 {
739 scratch_data.resize_system_vectors(n_components);
740
741 // So we have a multi-component DoFHandler here. That's more
742 // complicated. If the vector is real-valued, then we can just
743 // get everything at all quadrature points and copy them into
744 // the output array. In fact, we don't have to worry at all
745 // about the interpretation of the components.
746 if (dataset->is_complex_valued() == false)
747 {
748 dataset->get_function_values(
749 this_fe_patch_values,
751 real_part,
752 scratch_data.patch_values_system.solution_values);
753 for (unsigned int component = 0; component < n_components;
754 ++component)
755 for (unsigned int q = 0; q < n_q_points; ++q)
756 patch.data(offset + component, q) =
757 scratch_data.patch_values_system.solution_values[q](
758 component);
759
760 // Increment the counter for the actual data record.
761 offset += dataset->n_output_variables;
762 }
763 else
764 // The situation is more complicated if the input vector is
765 // complex-valued. The easiest approach would have been to
766 // just have all real and then all imaginary components.
767 // This would have been conceptually easy, but it has the
768 // annoying downside that if you have a vector-valued
769 // problem (say, [u v]) then the output order would have
770 // been [u_re, v_re, u_im, v_im]. That's tolerable, but not
771 // quite so nice because one typically thinks of real and
772 // imaginary parts as belonging together. We would really
773 // like the output order to be [u_re, u_im, v_re, v_im].
774 // That, too, would have been easy to implement because one
775 // just has to interleave real and imaginary parts.
776 //
777 // But that's also not what we want. That's because if one
778 // were, for example, to solve a complex-valued Stokes
779 // problem (e.g., computing eigenfunctions of the Stokes
780 // operator), then one has solution components
781 // [[u v] p] and the proper output order is
782 // [[u_re v_re] [u_im v_im] p_re p_im].
783 // In other words, the order in which we want to output
784 // data depends on the *interpretation* of components.
785 //
786 // Doing this requires a bit more code, and also needs to
787 // be in sync with what we do in
788 // DataOut_DoFData::get_dataset_names() and
789 // DataOut_DoFData::get_nonscalar_data_ranges().
790 {
791 // Given this description, first get the real parts of
792 // all components:
793 dataset->get_function_values(
794 this_fe_patch_values,
796 real_part,
797 scratch_data.patch_values_system.solution_values);
798
799 // Then we need to distribute them to the correct
800 // location. This requires knowledge of the interpretation
801 // of components as discussed above.
802 {
803 Assert(dataset->data_component_interpretation.size() ==
804 n_components,
806
807 unsigned int destination = offset;
808 for (unsigned int component = 0;
809 component < n_components;
810 /* component is updated below */)
811 {
812 switch (
813 dataset->data_component_interpretation[component])
814 {
817 {
818 // OK, a scalar component. Put all of the
819 // values into the current row
820 // ('destination'); then move 'component'
821 // forward by one (so we treat the next
822 // component) and 'destination' forward by
823 // two (because we're going to put the
824 // imaginary part of the current component
825 // into the next slot).
826 for (unsigned int q = 0; q < n_q_points;
827 ++q)
828 patch.data(destination, q) =
829 scratch_data.patch_values_system
830 .solution_values[q](component);
831
832 ++component;
833 destination += 2;
834
835 break;
836 }
837
840 {
841 // A vector component. Put the
842 // DoFHandlerType::space_dimension
843 // components into the next set of
844 // contiguous rows
845 // ('destination+c'); then move 'component'
846 // forward by spacedim (so we get to the
847 // next component after the current vector)
848 // and 'destination' forward by two*spacedim
849 // (because we're going to put the imaginary
850 // part of the vector into the subsequent
851 // spacedim slots).
852 const unsigned int size =
853 DoFHandlerType::space_dimension;
854 for (unsigned int c = 0; c < size; ++c)
855 for (unsigned int q = 0; q < n_q_points;
856 ++q)
857 patch.data(destination + c, q) =
858 scratch_data.patch_values_system
859 .solution_values[q](component + c);
860
861 component += size;
862 destination += 2 * size;
863
864 break;
865 }
866
869 {
870 // Same approach as for vectors above.
871 const unsigned int size =
872 DoFHandlerType::space_dimension *
873 DoFHandlerType::space_dimension;
874 for (unsigned int c = 0; c < size; ++c)
875 for (unsigned int q = 0; q < n_q_points;
876 ++q)
877 patch.data(destination + c, q) =
878 scratch_data.patch_values_system
879 .solution_values[q](component + c);
880
881 component += size;
882 destination += 2 * size;
883
884 break;
885 }
886
887 default:
888 Assert(false, ExcNotImplemented());
889 }
890 }
891 }
892
893 // And now we need to do the same thing again for the
894 // imaginary parts, starting at the top of the list of
895 // components/destinations again.
896 dataset->get_function_values(
897 this_fe_patch_values,
899 imaginary_part,
900 scratch_data.patch_values_system.solution_values);
901 {
902 unsigned int destination = offset;
903 for (unsigned int component = 0;
904 component < n_components;
905 /* component is updated below */)
906 {
907 switch (
908 dataset->data_component_interpretation[component])
909 {
912 {
913 // OK, a scalar component. Put all of the
914 // values into the row past the current one
915 // ('destination+1') since 'destination' is
916 // occupied by the real part.
917 for (unsigned int q = 0; q < n_q_points;
918 ++q)
919 patch.data(destination + 1, q) =
920 scratch_data.patch_values_system
921 .solution_values[q](component);
922
923 ++component;
924 destination += 2;
925
926 break;
927 }
928
931 {
932 // A vector component. Put the
933 // DoFHandlerType::space_dimension
934 // components into the set of contiguous
935 // rows that follow the real parts
936 // ('destination+spacedim+c').
937 const unsigned int size =
938 DoFHandlerType::space_dimension;
939 for (unsigned int c = 0; c < size; ++c)
940 for (unsigned int q = 0; q < n_q_points;
941 ++q)
942 patch.data(destination + size + c, q) =
943 scratch_data.patch_values_system
944 .solution_values[q](component + c);
945
946 component += size;
947 destination += 2 * size;
948
949 break;
950 }
951
954 {
955 // Same as for vectors.
956 const unsigned int size =
957 DoFHandlerType::space_dimension *
958 DoFHandlerType::space_dimension;
959 for (unsigned int c = 0; c < size; ++c)
960 for (unsigned int q = 0; q < n_q_points;
961 ++q)
962 patch.data(destination + size + c, q) =
963 scratch_data.patch_values_system
964 .solution_values[q](component + c);
965
966 component += size;
967 destination += 2 * size;
968
969 break;
970 }
971
972 default:
973 Assert(false, ExcNotImplemented());
974 }
975 }
976 }
977
978 // Increment the counter for the actual data record. We
979 // need to move it forward a number of positions equal to
980 // the number of components of this data set, times two
981 // because we dealt with a complex-valued input vector
982 offset += dataset->n_output_variables * 2;
983 }
984 }
985 }
986
987 // Also update the dataset_number index that we carry along with the
988 // for-loop over all data sets.
989 ++dataset_number;
990 }
991
992 // Then do the cell data. At least, we don't have to worry about
993 // complex-valued vectors/tensors since cell data is always scalar.
994 if (this->cell_data.size() != 0)
995 {
996 Assert(!cell_and_index->first->has_children(), ExcNotImplemented());
997
998 for (const auto &dataset : this->cell_data)
999 {
1000 // as above, first output the real part
1001 {
1002 const double value =
1003 dataset->get_cell_data_value(cell_and_index->second,
1005 ComponentExtractor::real_part);
1006 for (unsigned int q = 0; q < n_q_points; ++q)
1007 patch.data(offset, q) = value;
1008 }
1009
1010 // and if there is one, also output the imaginary part
1011 if (dataset->is_complex_valued() == true)
1012 {
1013 const double value = dataset->get_cell_data_value(
1014 cell_and_index->second,
1016 imaginary_part);
1017 for (unsigned int q = 0; q < n_q_points; ++q)
1018 patch.data(offset + 1, q) = value;
1019 }
1020
1021 offset += (dataset->is_complex_valued() ? 2 : 1);
1022 }
1023 }
1024 }
1025
1026
1027 for (const unsigned int f : cell_and_index->first->face_indices())
1028 {
1029 // let's look up whether the neighbor behind that face is noted in the
1030 // table of cells which we treat. this can only happen if the neighbor
1031 // exists, and is on the same level as this cell, but it may also happen
1032 // that the neighbor is not a member of the range of cells over which we
1033 // loop, in which case the respective entry in the
1034 // cell_to_patch_index_map will have the value no_neighbor. (note that
1035 // since we allocated only as much space in this array as the maximum
1036 // index of the cells we loop over, not every neighbor may have its
1037 // space in it, so we have to assume that it is extended by values
1038 // no_neighbor)
1039 if (cell_and_index->first->at_boundary(f) ||
1040 (cell_and_index->first->neighbor(f)->level() !=
1041 cell_and_index->first->level()))
1042 {
1043 patch.neighbors[f] = numbers::invalid_unsigned_int;
1044 continue;
1045 }
1046
1047 const cell_iterator neighbor = cell_and_index->first->neighbor(f);
1048 Assert(static_cast<unsigned int>(neighbor->level()) <
1049 scratch_data.cell_to_patch_index_map->size(),
1051 if ((static_cast<unsigned int>(neighbor->index()) >=
1052 (*scratch_data.cell_to_patch_index_map)[neighbor->level()].size()) ||
1053 ((*scratch_data.cell_to_patch_index_map)[neighbor->level()]
1054 [neighbor->index()] ==
1056 {
1057 patch.neighbors[f] = numbers::invalid_unsigned_int;
1058 continue;
1059 }
1060
1061 // now, there is a neighbor, so get its patch number and set it for the
1062 // neighbor index
1063 patch.neighbors[f] =
1064 (*scratch_data
1065 .cell_to_patch_index_map)[neighbor->level()][neighbor->index()];
1066 }
1067
1068 const unsigned int patch_idx =
1069 (*scratch_data.cell_to_patch_index_map)[cell_and_index->first->level()]
1070 [cell_and_index->first->index()];
1071 // did we mess up the indices?
1072 Assert(patch_idx < this->patches.size(), ExcInternalError());
1073 patch.patch_index = patch_idx;
1074
1075 // Put the patch into the patches vector. instead of copying the data,
1076 // simply swap the contents to avoid the penalty of writing into another
1077 // processor's memory
1078 this->patches[patch_idx].swap(patch);
1079}
1080
1081
1082
1083template <int dim, typename DoFHandlerType>
1084void
1085DataOut<dim, DoFHandlerType>::build_patches(const unsigned int n_subdivisions)
1086{
1087 AssertDimension(this->triangulation->get_reference_cells().size(), 1);
1088
1089 build_patches(
1090 this->triangulation->get_reference_cells()[0]
1091 .template get_default_linear_mapping<DoFHandlerType::dimension,
1092 DoFHandlerType::space_dimension>(),
1093 n_subdivisions,
1094 no_curved_cells);
1095}
1096
1097
1098
1099template <int dim, typename DoFHandlerType>
1100void
1103 & mapping,
1104 const unsigned int n_subdivisions_,
1105 const CurvedCellRegion curved_region)
1106{
1107 hp::MappingCollection<DoFHandlerType::dimension,
1108 DoFHandlerType::space_dimension>
1109 mapping_collection(mapping);
1110
1111 build_patches(mapping_collection, n_subdivisions_, curved_region);
1112}
1113
1114
1115
1116template <int dim, typename DoFHandlerType>
1117void
1119 const hp::MappingCollection<DoFHandlerType::dimension,
1120 DoFHandlerType::space_dimension> &mapping,
1121 const unsigned int n_subdivisions_,
1122 const CurvedCellRegion curved_region)
1123{
1124 // Check consistency of redundant template parameter
1125 Assert(dim == DoFHandlerType::dimension,
1126 ExcDimensionMismatch(dim, DoFHandlerType::dimension));
1127
1128 Assert(this->triangulation != nullptr,
1130
1131 const unsigned int n_subdivisions =
1132 (n_subdivisions_ != 0) ? n_subdivisions_ : this->default_subdivisions;
1133 Assert(n_subdivisions >= 1,
1135 n_subdivisions));
1136
1137 this->validate_dataset_names();
1138
1139 // First count the cells we want to create patches of. Also fill the object
1140 // that maps the cell indices to the patch numbers, as this will be needed
1141 // for generation of neighborship information.
1142 // Note, there is a confusing mess of different indices here at play:
1143 // - patch_index: the index of a patch in all_cells
1144 // - cell->index: only unique on each level, used in cell_to_patch_index_map
1145 // - active_index: index for a cell when counting from begin_active() using
1146 // ++cell (identical to cell->active_cell_index())
1147 // - cell_index: unique index of a cell counted using
1148 // next_cell_function() starting from first_cell_function()
1149 //
1150 // It turns out that we create one patch for each selected cell, so
1151 // patch_index==cell_index.
1152 //
1153 // Now construct the map such that
1154 // cell_to_patch_index_map[cell->level][cell->index] = patch_index
1155 std::vector<std::vector<unsigned int>> cell_to_patch_index_map;
1156 cell_to_patch_index_map.resize(this->triangulation->n_levels());
1157 for (unsigned int l = 0; l < this->triangulation->n_levels(); ++l)
1158 {
1159 // max_index is the largest cell->index on level l
1160 unsigned int max_index = 0;
1161 for (cell_iterator cell = first_cell_function(*this->triangulation);
1162 cell != this->triangulation->end();
1163 cell = next_cell_function(*this->triangulation, cell))
1164 if (static_cast<unsigned int>(cell->level()) == l)
1165 max_index =
1166 std::max(max_index, static_cast<unsigned int>(cell->index()));
1167
1168 cell_to_patch_index_map[l].resize(
1169 max_index + 1,
1171 DoFHandlerType::dimension,
1172 DoFHandlerType::space_dimension>::no_neighbor);
1173 }
1174
1175 // will be all_cells[patch_index] = pair(cell, active_index)
1176 std::vector<std::pair<cell_iterator, unsigned int>> all_cells;
1177 {
1178 // important: we need to compute the active_index of the cell in the range
1179 // 0..n_active_cells() because this is where we need to look up cell
1180 // data from (cell data vectors do not have the length distance computed by
1181 // first_cell_function/next_cell_function because this might skip
1182 // some values (FilteredIterator).
1183 auto active_cell = this->triangulation->begin_active();
1184 unsigned int active_index = 0;
1185 cell_iterator cell = first_cell_function(*this->triangulation);
1186 for (; cell != this->triangulation->end();
1187 cell = next_cell_function(*this->triangulation, cell))
1188 {
1189 // move forward until active_cell points at the cell (cell) we are
1190 // looking at to compute the current active_index
1191 while (active_cell != this->triangulation->end() && cell->is_active() &&
1192 decltype(active_cell)(cell) != active_cell)
1193 {
1194 ++active_cell;
1195 ++active_index;
1196 }
1197
1198 Assert(static_cast<unsigned int>(cell->level()) <
1199 cell_to_patch_index_map.size(),
1201 Assert(static_cast<unsigned int>(cell->index()) <
1202 cell_to_patch_index_map[cell->level()].size(),
1204 Assert(active_index < this->triangulation->n_active_cells(),
1206 cell_to_patch_index_map[cell->level()][cell->index()] =
1207 all_cells.size();
1208
1209 all_cells.emplace_back(cell, active_index);
1210 }
1211 }
1212
1213 this->patches.clear();
1214 this->patches.resize(all_cells.size());
1215
1216 // Now create a default object for the WorkStream object to work with. The
1217 // first step is to count how many output data sets there will be. This is,
1218 // in principle, just the number of components of each data set, but we
1219 // need to allocate two entries per component if there are
1220 // complex-valued input data (unless we use a postprocessor on this
1221 // output -- all postprocessor outputs are real-valued)
1222 unsigned int n_datasets = 0;
1223 for (unsigned int i = 0; i < this->cell_data.size(); ++i)
1224 n_datasets += (this->cell_data[i]->is_complex_valued() &&
1225 (this->cell_data[i]->postprocessor == nullptr) ?
1226 2 :
1227 1);
1228 for (unsigned int i = 0; i < this->dof_data.size(); ++i)
1229 n_datasets += (this->dof_data[i]->n_output_variables *
1230 (this->dof_data[i]->is_complex_valued() &&
1231 (this->dof_data[i]->postprocessor == nullptr) ?
1232 2 :
1233 1));
1234
1235 std::vector<unsigned int> n_postprocessor_outputs(this->dof_data.size());
1236 for (unsigned int dataset = 0; dataset < this->dof_data.size(); ++dataset)
1237 if (this->dof_data[dataset]->postprocessor)
1238 n_postprocessor_outputs[dataset] =
1239 this->dof_data[dataset]->n_output_variables;
1240 else
1241 n_postprocessor_outputs[dataset] = 0;
1242
1243 const CurvedCellRegion curved_cell_region =
1244 (n_subdivisions < 2 ? no_curved_cells : curved_region);
1245
1246 UpdateFlags update_flags = update_values;
1247 if (curved_cell_region != no_curved_cells)
1248 update_flags |= update_quadrature_points;
1249
1250 for (unsigned int i = 0; i < this->dof_data.size(); ++i)
1251 if (this->dof_data[i]->postprocessor)
1252 update_flags |=
1253 this->dof_data[i]->postprocessor->get_needed_update_flags();
1254 // perhaps update_normal_vectors is present, which would only be useful on
1255 // faces, but we may not use it here.
1256 Assert(
1257 !(update_flags & update_normal_vectors),
1258 ExcMessage(
1259 "The update of normal vectors may not be requested for evaluation of "
1260 "data on cells via DataPostprocessor."));
1261
1262 internal::DataOutImplementation::ParallelData<DoFHandlerType::dimension,
1263 DoFHandlerType::space_dimension>
1264 thread_data(n_datasets,
1265 n_subdivisions,
1266 n_postprocessor_outputs,
1267 mapping,
1268 this->get_fes(),
1269 update_flags,
1270 cell_to_patch_index_map);
1271
1272 auto worker = [this, n_subdivisions, curved_cell_region](
1273 const std::pair<cell_iterator, unsigned int> *cell_and_index,
1275 DoFHandlerType::dimension,
1276 DoFHandlerType::space_dimension> &scratch_data,
1277 // this function doesn't actually need a copy data object --
1278 // it just writes everything right into the output array
1279 int) {
1280 this->build_one_patch(cell_and_index,
1281 scratch_data,
1282 n_subdivisions,
1283 curved_cell_region);
1284 };
1285
1286 // now build the patches in parallel
1287 if (all_cells.size() > 0)
1288 WorkStream::run(all_cells.data(),
1289 all_cells.data() + all_cells.size(),
1290 worker,
1291 // no copy-local-to-global function needed here
1292 std::function<void(const int)>(),
1293 thread_data,
1294 /* dummy CopyData object = */ 0,
1295 // experimenting shows that we can make things run a bit
1296 // faster if we increase the number of cells we work on
1297 // per item (i.e., WorkStream's chunk_size argument,
1298 // about 10% improvement) and the items in flight at any
1299 // given time (another 5% on the testcase discussed in
1300 // @ref workstream_paper, on 32 cores) and if
1302 64);
1303}
1304
1305
1306
1307template <int dim, typename DoFHandlerType>
1308void
1310 const std::function<cell_iterator(const Triangulation<dim, spacedim> &)>
1311 & first_cell,
1312 const std::function<cell_iterator(const Triangulation<dim, spacedim> &,
1313 const cell_iterator &)> &next_cell)
1314{
1315 first_cell_function = first_cell;
1316 next_cell_function = next_cell;
1317}
1318
1319
1320
1321template <int dim, typename DoFHandlerType>
1322void
1324 const FilteredIterator<cell_iterator> &filtered_iterator)
1325{
1326 const auto first_cell =
1327 [filtered_iterator](const Triangulation<dim, spacedim> &triangulation) {
1328 // Create a copy of the filtered iterator so that we can
1329 // call a non-const function -- though we are really only
1330 // interested in the return value of that function, not the
1331 // state of the object
1332 FilteredIterator<cell_iterator> x = filtered_iterator;
1333 return x.set_to_next_positive(triangulation.begin());
1334 };
1335
1336
1337 const auto next_cell =
1338 [filtered_iterator](const Triangulation<dim, spacedim> &,
1339 const cell_iterator &cell) {
1340 // Create a copy of the filtered iterator so that we can
1341 // call a non-const function -- though we are really only
1342 // interested in the return value of that function, not the
1343 // state of the object
1344 FilteredIterator<cell_iterator> x = filtered_iterator;
1345
1346 // Set the iterator to 'cell'. Since 'cell' must satisfy the
1347 // predicate (that's how it was created), set_to_next_positive
1348 // simply sets the iterator to 'cell'.
1349 x.set_to_next_positive(cell);
1350
1351 // Advance by one:
1352 ++x;
1353
1354 return x;
1355 };
1356
1357 set_cell_selection(first_cell, next_cell);
1358}
1359
1360
1361
1362template <int dim, typename DoFHandlerType>
1363const std::pair<typename DataOut<dim, DoFHandlerType>::FirstCellFunctionType,
1366{
1367 return std::make_pair(first_cell_function, next_cell_function);
1368}
1369
1370
1371
1372template <int dim, typename DoFHandlerType>
1375{
1376 return this->triangulation->begin_active();
1377}
1378
1379
1380
1381template <int dim, typename DoFHandlerType>
1385{
1386 // convert the iterator to an active_iterator and advance this to the next
1387 // active cell
1388 typename Triangulation<DoFHandlerType::dimension,
1389 DoFHandlerType::space_dimension>::active_cell_iterator
1390 active_cell = cell;
1391 ++active_cell;
1392 return active_cell;
1393}
1394
1395
1396
1397template <int dim, typename DoFHandlerType>
1400{
1401 typename DataOut<dim, DoFHandlerType>::cell_iterator cell = first_cell();
1402
1403 // skip cells if the current one has no children (is active) and is a ghost
1404 // or artificial cell
1405 while ((cell != this->triangulation->end()) && cell->is_active() &&
1406 !cell->is_locally_owned())
1407 cell = next_cell(cell);
1408
1409 return cell;
1410}
1411
1412
1413
1414template <int dim, typename DoFHandlerType>
1417 const typename DataOut<dim, DoFHandlerType>::cell_iterator &old_cell)
1418{
1420 next_cell(old_cell);
1421 while ((cell != this->triangulation->end()) && cell->is_active() &&
1422 !cell->is_locally_owned())
1423 cell = next_cell(cell);
1424 return cell;
1425}
1426
1427
1428// explicit instantiations
1429#include "data_out.inst"
1430
CurvedCellRegion
Definition: data_out.h:196
virtual cell_iterator first_cell()
Definition: data_out.cc:1374
typename DataOut_DoFData< DoFHandlerType, DoFHandlerType::dimension, DoFHandlerType::space_dimension >::cell_iterator cell_iterator
Definition: data_out.h:165
void set_cell_selection(const std::function< cell_iterator(const Triangulation< dim, spacedim > &)> &first_cell, const std::function< cell_iterator(const Triangulation< dim, spacedim > &, const cell_iterator &)> &next_cell)
Definition: data_out.cc:1309
virtual void build_patches(const unsigned int n_subdivisions=0)
Definition: data_out.cc:1085
const std::pair< FirstCellFunctionType, NextCellFunctionType > get_cell_selection() const
Definition: data_out.cc:1365
virtual cell_iterator first_locally_owned_cell()
Definition: data_out.cc:1399
typename std::function< cell_iterator(const Triangulation< dim, spacedim > &, const cell_iterator &)> NextCellFunctionType
Definition: data_out.h:180
DataOut()
Definition: data_out.cc:69
void build_one_patch(const std::pair< cell_iterator, unsigned int > *cell_and_index, internal::DataOutImplementation::ParallelData< DoFHandlerType::dimension, DoFHandlerType::space_dimension > &scratch_data, const unsigned int n_subdivisions, const CurvedCellRegion curved_cell_region)
Definition: data_out.cc:88
virtual cell_iterator next_locally_owned_cell(const cell_iterator &cell)
Definition: data_out.cc:1416
virtual cell_iterator next_cell(const cell_iterator &cell)
Definition: data_out.cc:1383
virtual UpdateFlags get_needed_update_flags() const =0
virtual void evaluate_vector_field(const DataPostprocessorInputs::Vector< dim > &input_data, std::vector< Vector< double > > &computed_quantities) const
virtual void evaluate_scalar_field(const DataPostprocessorInputs::Scalar< dim > &input_data, std::vector< Vector< double > > &computed_quantities) const
const std::vector< Point< spacedim > > & get_quadrature_points() const
const Mapping< dim, spacedim > & get_mapping() const
const unsigned int n_quadrature_points
Definition: fe_values.h:2432
FilteredIterator & set_to_next_positive(const BaseIterator &bi)
Abstract base class for mapping classes.
Definition: mapping.h:304
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const =0
static unsigned int n_threads()
Definition: tensor.h:472
Definition: vector.h:110
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
UpdateFlags
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
static ::ExceptionBase & ExcInvalidNumberOfSubdivisions(int arg1)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1465
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
static ::ExceptionBase & ExcInternalError()
unsigned int n_subdivisions
static ::ExceptionBase & ExcNoTriangulationSelected()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
Definition: mapping.cc:240
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void run(const std::vector< std::vector< Iterator > > &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: work_stream.h:1252
static const unsigned int invalid_unsigned_int
Definition: types.h:196
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
ParallelData(const unsigned int n_datasets, const unsigned int n_subdivisions, const std::vector< unsigned int > &n_postprocessor_outputs, const ::hp::MappingCollection< dim, spacedim > &mapping, const std::vector< std::shared_ptr<::hp::FECollection< dim, spacedim > > > &finite_elements, const UpdateFlags update_flags, const std::vector< std::vector< unsigned int > > &cell_to_patch_index_map)
Definition: data_out.cc:44