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