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