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