Reference documentation for deal.II version GIT 0b65fff18a 2023-09-27 19:30:02+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\}}\)
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>
23 #include <deal.II/fe/fe_values.h>
24 #include <deal.II/fe/mapping.h>
25 
26 #include <deal.II/grid/tria.h>
28 
29 #include <deal.II/hp/fe_values.h>
30 
32 
33 #include <sstream>
34 
36 
37 
38 namespace 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 
67 template <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 
94 template <int dim, int spacedim>
95 void
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,
353  ComponentExtractor::real_part,
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);
367  scratch_data.patch_values_system
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,
436  ComponentExtractor::imaginary_part,
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,
497  ComponentExtractor::real_part,
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,
524  ComponentExtractor::imaginary_part,
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,
552  ComponentExtractor::real_part,
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,
575  ComponentExtractor::imaginary_part,
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,
603  ComponentExtractor::real_part,
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,
626  ComponentExtractor::imaginary_part,
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,
711  imaginary_part,
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,
787  ExcInternalError());
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,
879  imaginary_part,
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,
983  ComponentExtractor::real_part);
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,
994  imaginary_part);
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(),
1028  ExcInternalError());
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 
1061 template <int dim, int spacedim>
1062 void
1063 DataOut<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 
1075 template <int dim, int spacedim>
1076 void
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 
1088 template <int dim, int spacedim>
1089 void
1091  const hp::MappingCollection<dim, spacedim> &mapping,
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(),
1167  ExcInternalError());
1168  Assert(static_cast<unsigned int>(cell->index()) <
1169  cell_to_patch_index_map[cell->level()].size(),
1170  ExcInternalError());
1171  Assert(active_index < this->triangulation->n_active_cells(),
1172  ExcInternalError());
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 
1272 template <int dim, int spacedim>
1273 void
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 
1286 template <int dim, int spacedim>
1287 void
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 
1327 template <int dim, int spacedim>
1328 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
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
std::pair< FirstCellFunctionType, NextCellFunctionType > get_cell_selection() const
Definition: data_out.cc:1330
virtual void evaluate_scalar_field(const DataPostprocessorInputs::Scalar< dim > &input_data, std::vector< Vector< double >> &computed_quantities) const
virtual void evaluate_vector_field(const DataPostprocessorInputs::Vector< dim > &input_data, std::vector< Vector< double >> &computed_quantities) const
virtual UpdateFlags get_needed_update_flags() const =0
const unsigned int n_quadrature_points
const Mapping< dim, spacedim > & get_mapping() const
const FiniteElement< dim, spacedim > & get_fe() const
const std::vector< Point< spacedim > > & get_quadrature_points() 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:110
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
Point< 3 > vertices[4]
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcInvalidNumberOfSubdivisions(int arg1)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNoTriangulationSelected()
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1789
static ::ExceptionBase & ExcMessage(std::string arg1)
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.
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:1270
void copy(const T *begin, const T *end, U *dest)
static const unsigned int invalid_unsigned_int
Definition: types.h:213
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
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)
const FEValuesBase< dim, spacedim > & get_present_fe_values(const unsigned int dataset) const
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