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