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