Reference documentation for deal.II version Git 0491297983 2019-09-23 09:31:59 +0200
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
data_out.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2019 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 
16 #include <deal.II/base/work_stream.h>
17 
18 #include <deal.II/dofs/dof_accessor.h>
19 #include <deal.II/dofs/dof_handler.h>
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>
27 #include <deal.II/grid/tria_iterator.h>
28 
29 #include <deal.II/hp/dof_handler.h>
30 #include <deal.II/hp/fe_values.h>
31 
32 #include <deal.II/numerics/data_out.h>
33 
34 #include <sstream>
35 
36 DEAL_II_NAMESPACE_OPEN
37 
38 
39 namespace internal
40 {
41  namespace DataOutImplementation
42  {
43  template <int dim, int spacedim>
44  ParallelData<dim, spacedim>::ParallelData(
45  const unsigned int n_datasets,
46  const unsigned int n_subdivisions,
47  const std::vector<unsigned int> &n_postprocessor_outputs,
48  const Mapping<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, typename DoFHandlerType>
69 void
71  const std::pair<cell_iterator, unsigned int> *cell_and_index,
72  internal::DataOutImplementation::ParallelData<DoFHandlerType::dimension,
73  DoFHandlerType::space_dimension>
74  & scratch_data,
75  const unsigned int n_subdivisions,
76  const CurvedCellRegion curved_cell_region)
77 {
78  // first create the output object that we will write into
79  ::DataOutBase::Patch<DoFHandlerType::dimension,
80  DoFHandlerType::space_dimension>
81  patch;
82  patch.n_subdivisions = n_subdivisions;
83 
84  // set the vertices of the patch. if the mapping does not preserve locations
85  // (e.g. MappingQEulerian), we need to compute the offset of the vertex for
86  // the graphical output. Otherwise, we can just use the vertex info.
87  for (unsigned int vertex = 0;
88  vertex < GeometryInfo<DoFHandlerType::dimension>::vertices_per_cell;
89  ++vertex)
90  if (scratch_data.mapping_collection[0].preserves_vertex_locations())
91  patch.vertices[vertex] = cell_and_index->first->vertex(vertex);
92  else
93  patch.vertices[vertex] =
94  scratch_data.mapping_collection[0].transform_unit_to_real_cell(
95  cell_and_index->first,
97 
98  // create DoFHandlerType::active_cell_iterator and initialize FEValues
99  scratch_data.reinit_all_fe_values(this->dof_data, cell_and_index->first);
100 
102  &fe_patch_values = scratch_data.get_present_fe_values(0);
103 
104  const unsigned int n_q_points = fe_patch_values.n_quadrature_points;
105 
106  // depending on the requested output of curved cells, if necessary
107  // append the quadrature points to the last rows of the patch.data
108  // member. This is the case if we want to produce curved cells at the
109  // boundary and this cell actually is at the boundary, or else if we
110  // want to produce curved cells everywhere
111  //
112  // note: a cell is *always* at the boundary if dim<spacedim
113  if (curved_cell_region == curved_inner_cells ||
114  (curved_cell_region == curved_boundary &&
115  (cell_and_index->first->at_boundary() ||
116  (DoFHandlerType::dimension != DoFHandlerType::space_dimension))))
117  {
118  Assert(patch.space_dim == DoFHandlerType::space_dimension,
119  ExcInternalError());
120 
121  // set the flag indicating that for this cell the points are
122  // explicitly given
123  patch.points_are_available = true;
124 
125  // then resize the patch.data member in order to have enough memory for
126  // the quadrature points as well, and copy the quadrature points there
127  const std::vector<Point<DoFHandlerType::space_dimension>> &q_points =
128  fe_patch_values.get_quadrature_points();
129 
130  patch.data.reinit(scratch_data.n_datasets +
131  DoFHandlerType::space_dimension,
132  n_q_points);
133  for (unsigned int i = 0; i < DoFHandlerType::space_dimension; ++i)
134  for (unsigned int q = 0; q < n_q_points; ++q)
135  patch.data(patch.data.size(0) - DoFHandlerType::space_dimension + i,
136  q) = q_points[q][i];
137  }
138  else
139  {
140  patch.data.reinit(scratch_data.n_datasets, n_q_points);
141  patch.points_are_available = false;
142  }
143 
144 
145  if (scratch_data.n_datasets > 0)
146  {
147  // counter for data records
148  unsigned int offset = 0;
149 
150  // first fill dof_data
151  for (unsigned int dataset = 0; dataset < this->dof_data.size(); ++dataset)
152  {
153  const FEValuesBase<DoFHandlerType::dimension,
154  DoFHandlerType::space_dimension>
155  &this_fe_patch_values = scratch_data.get_present_fe_values(dataset);
156  const unsigned int n_components =
157  this_fe_patch_values.get_fe().n_components();
158 
160  *postprocessor = this->dof_data[dataset]->postprocessor;
161 
162  if (postprocessor != nullptr)
163  {
164  // We have to postprocess the data, so determine, which fields
165  // have to be updated
166  const UpdateFlags update_flags =
167  postprocessor->get_needed_update_flags();
168 
169  if ((n_components == 1) &&
170  (this->dof_data[dataset]->is_complex_valued() == false))
171  {
172  // At each point there is only one component of value,
173  // gradient etc. Based on the 'if' statement above, we
174  // know that the solution is scalar and real-valued, so
175  // we do not need to worry about getting any imaginary
176  // components to the postprocessor, and we can safely
177  // call the function that evaluates a scalar field
178  if (update_flags & update_values)
179  this->dof_data[dataset]->get_function_values(
180  this_fe_patch_values,
181  internal::DataOutImplementation::ComponentExtractor::
182  real_part,
183  scratch_data.patch_values_scalar.solution_values);
184 
185  if (update_flags & update_gradients)
186  this->dof_data[dataset]->get_function_gradients(
187  this_fe_patch_values,
188  internal::DataOutImplementation::ComponentExtractor::
189  real_part,
190  scratch_data.patch_values_scalar.solution_gradients);
191 
192  if (update_flags & update_hessians)
193  this->dof_data[dataset]->get_function_hessians(
194  this_fe_patch_values,
195  internal::DataOutImplementation::ComponentExtractor::
196  real_part,
197  scratch_data.patch_values_scalar.solution_hessians);
198 
199  // Also fill some of the other fields postprocessors may
200  // want to access.
201  if (update_flags & update_quadrature_points)
202  scratch_data.patch_values_scalar.evaluation_points =
203  this_fe_patch_values.get_quadrature_points();
204 
205  const typename DoFHandlerType::active_cell_iterator dh_cell(
206  &cell_and_index->first->get_triangulation(),
207  cell_and_index->first->level(),
208  cell_and_index->first->index(),
209  this->dof_data[dataset]->dof_handler);
210  scratch_data.patch_values_scalar
211  .template set_cell<DoFHandlerType>(dh_cell);
212 
213  // Finally call the postprocessor's function that
214  // deals with scalar inputs.
215  postprocessor->evaluate_scalar_field(
216  scratch_data.patch_values_scalar,
217  scratch_data.postprocessed_values[dataset]);
218  }
219  else
220  {
221  scratch_data.resize_system_vectors(
222  n_components *
223  (this->dof_data[dataset]->is_complex_valued() ? 2 : 1));
224 
225  // At each point we now have to evaluate a vector valued
226  // function and its derivatives. It may be that the solution
227  // is scalar and complex-valued, but we treat this as a vector
228  // field with two components.
229 
230  // At this point, we need to ask how we fill the fields that
231  // we want to pass on to the postprocessor. If the field in
232  // question is real-valued, we'll just extract the (only)
233  // real component from the solution fields
234  if (this->dof_data[dataset]->is_complex_valued() == false)
235  {
236  if (update_flags & update_values)
237  this->dof_data[dataset]->get_function_values(
238  this_fe_patch_values,
239  internal::DataOutImplementation::ComponentExtractor::
240  real_part,
241  scratch_data.patch_values_system.solution_values);
242 
243  if (update_flags & update_gradients)
244  this->dof_data[dataset]->get_function_gradients(
245  this_fe_patch_values,
246  internal::DataOutImplementation::ComponentExtractor::
247  real_part,
248  scratch_data.patch_values_system.solution_gradients);
249 
250  if (update_flags & update_hessians)
251  this->dof_data[dataset]->get_function_hessians(
252  this_fe_patch_values,
253  internal::DataOutImplementation::ComponentExtractor::
254  real_part,
255  scratch_data.patch_values_system.solution_hessians);
256  }
257  else
258  {
259  // The solution is complex-valued. We don't currently
260  // know how to handle this in the most general case,
261  // but we can deal with it as long as there is only a
262  // scalar solution since then we can just collate the two
263  // components of the scalar solution into one vector field
264  if (n_components == 1)
265  {
266  // First get the real component of the scalar solution
267  // and copy the data into the
268  // scratch_data.patch_values_system output fields
269  if (update_flags & update_values)
270  {
271  this->dof_data[dataset]->get_function_values(
272  this_fe_patch_values,
274  ComponentExtractor::real_part,
275  scratch_data.patch_values_scalar
276  .solution_values);
277 
278  for (unsigned int i = 0;
279  i < scratch_data.patch_values_scalar
280  .solution_values.size();
281  ++i)
282  {
284  scratch_data.patch_values_system
285  .solution_values[i]
286  .size(),
287  2);
288  scratch_data.patch_values_system
289  .solution_values[i][0] =
290  scratch_data.patch_values_scalar
291  .solution_values[i];
292  }
293  }
294 
295  if (update_flags & update_gradients)
296  {
297  this->dof_data[dataset]->get_function_gradients(
298  this_fe_patch_values,
300  ComponentExtractor::real_part,
301  scratch_data.patch_values_scalar
302  .solution_gradients);
303 
304  for (unsigned int i = 0;
305  i < scratch_data.patch_values_scalar
306  .solution_gradients.size();
307  ++i)
308  {
310  scratch_data.patch_values_system
311  .solution_values[i]
312  .size(),
313  2);
314  scratch_data.patch_values_system
315  .solution_gradients[i][0] =
316  scratch_data.patch_values_scalar
317  .solution_gradients[i];
318  }
319  }
320 
321  if (update_flags & update_hessians)
322  {
323  this->dof_data[dataset]->get_function_hessians(
324  this_fe_patch_values,
326  ComponentExtractor::real_part,
327  scratch_data.patch_values_scalar
328  .solution_hessians);
329 
330  for (unsigned int i = 0;
331  i < scratch_data.patch_values_scalar
332  .solution_hessians.size();
333  ++i)
334  {
336  scratch_data.patch_values_system
337  .solution_hessians[i]
338  .size(),
339  2);
340  scratch_data.patch_values_system
341  .solution_hessians[i][0] =
342  scratch_data.patch_values_scalar
343  .solution_hessians[i];
344  }
345  }
346 
347  // Now we also have to get the imaginary
348  // component of the scalar solution
349  // and copy the data into the
350  // scratch_data.patch_values_system output fields
351  // that follow the real one
352  if (update_flags & update_values)
353  {
354  this->dof_data[dataset]->get_function_values(
355  this_fe_patch_values,
357  ComponentExtractor::imaginary_part,
358  scratch_data.patch_values_scalar
359  .solution_values);
360 
361  for (unsigned int i = 0;
362  i < scratch_data.patch_values_scalar
363  .solution_values.size();
364  ++i)
365  {
367  scratch_data.patch_values_system
368  .solution_values[i]
369  .size(),
370  2);
371  scratch_data.patch_values_system
372  .solution_values[i][1] =
373  scratch_data.patch_values_scalar
374  .solution_values[i];
375  }
376  }
377 
378  if (update_flags & update_gradients)
379  {
380  this->dof_data[dataset]->get_function_gradients(
381  this_fe_patch_values,
383  ComponentExtractor::imaginary_part,
384  scratch_data.patch_values_scalar
385  .solution_gradients);
386 
387  for (unsigned int i = 0;
388  i < scratch_data.patch_values_scalar
389  .solution_gradients.size();
390  ++i)
391  {
393  scratch_data.patch_values_system
394  .solution_values[i]
395  .size(),
396  2);
397  scratch_data.patch_values_system
398  .solution_gradients[i][1] =
399  scratch_data.patch_values_scalar
400  .solution_gradients[i];
401  }
402  }
403 
404  if (update_flags & update_hessians)
405  {
406  this->dof_data[dataset]->get_function_hessians(
407  this_fe_patch_values,
409  ComponentExtractor::imaginary_part,
410  scratch_data.patch_values_scalar
411  .solution_hessians);
412 
413  for (unsigned int i = 0;
414  i < scratch_data.patch_values_scalar
415  .solution_hessians.size();
416  ++i)
417  {
419  scratch_data.patch_values_system
420  .solution_hessians[i]
421  .size(),
422  2);
423  scratch_data.patch_values_system
424  .solution_hessians[i][1] =
425  scratch_data.patch_values_scalar
426  .solution_hessians[i];
427  }
428  }
429  }
430  else
431  {
432  // This is the vector-valued, complex-valued case.
433  Assert(false, ExcNotImplemented());
434  }
435  }
436 
437  // Now set other fields we may need
438  if (update_flags & update_quadrature_points)
439  scratch_data.patch_values_system.evaluation_points =
440  this_fe_patch_values.get_quadrature_points();
441 
442  const typename DoFHandlerType::active_cell_iterator dh_cell(
443  &cell_and_index->first->get_triangulation(),
444  cell_and_index->first->level(),
445  cell_and_index->first->index(),
446  this->dof_data[dataset]->dof_handler);
447  scratch_data.patch_values_system
448  .template set_cell<DoFHandlerType>(dh_cell);
449 
450  postprocessor->evaluate_vector_field(
451  scratch_data.patch_values_system,
452  scratch_data.postprocessed_values[dataset]);
453  }
454 
455  // Now we need to copy the result of the postprocessor to
456  // the Patch object where it can then be further processes
457  // by the functions in DataOutBase
458  for (unsigned int q = 0; q < n_q_points; ++q)
459  for (unsigned int component = 0;
460  component < this->dof_data[dataset]->n_output_variables;
461  ++component)
462  patch.data(offset + component, q) =
463  scratch_data.postprocessed_values[dataset][q](component);
464  }
465  else
466  {
467  // use the given data vector directly, without a postprocessor.
468  // again, we treat single component functions separately for
469  // efficiency reasons.
470  if (n_components == 1)
471  {
472  // first output the real part of the solution vector
473  this->dof_data[dataset]->get_function_values(
474  this_fe_patch_values,
475  internal::DataOutImplementation::ComponentExtractor::
476  real_part,
477  scratch_data.patch_values_scalar.solution_values);
478  for (unsigned int q = 0; q < n_q_points; ++q)
479  patch.data(offset, q) =
480  scratch_data.patch_values_scalar.solution_values[q];
481 
482  // and if there is one, also output the imaginary part
483  if (this->dof_data[dataset]->is_complex_valued() == true)
484  {
485  this->dof_data[dataset]->get_function_values(
486  this_fe_patch_values,
487  internal::DataOutImplementation::ComponentExtractor::
488  imaginary_part,
489  scratch_data.patch_values_scalar.solution_values);
490  for (unsigned int q = 0; q < n_q_points; ++q)
491  patch.data(offset + 1, q) =
492  scratch_data.patch_values_scalar.solution_values[q];
493  }
494  }
495  else
496  {
497  scratch_data.resize_system_vectors(n_components);
498 
499  // same as above: first the real part
500  const unsigned int stride =
501  (this->dof_data[dataset]->is_complex_valued() ? 2 : 1);
502  this->dof_data[dataset]->get_function_values(
503  this_fe_patch_values,
504  internal::DataOutImplementation::ComponentExtractor::
505  real_part,
506  scratch_data.patch_values_system.solution_values);
507  for (unsigned int component = 0; component < n_components;
508  ++component)
509  for (unsigned int q = 0; q < n_q_points; ++q)
510  patch.data(offset + component * stride, q) =
511  scratch_data.patch_values_system.solution_values[q](
512  component);
513 
514  // and if there is one, also output the imaginary part
515  if (this->dof_data[dataset]->is_complex_valued() == true)
516  {
517  this->dof_data[dataset]->get_function_values(
518  this_fe_patch_values,
519  internal::DataOutImplementation::ComponentExtractor::
520  imaginary_part,
521  scratch_data.patch_values_system.solution_values);
522  for (unsigned int component = 0; component < n_components;
523  ++component)
524  for (unsigned int q = 0; q < n_q_points; ++q)
525  patch.data(offset + component * stride + 1, q) =
526  scratch_data.patch_values_system.solution_values[q](
527  component);
528  }
529  }
530  }
531 
532  // Increment the counter for the actual data record. We need to
533  // move it forward a number of positions equal to the number
534  // of components of this data set; if the input consisted
535  // of a complex-valued quantity and if it is not further
536  // processed by a postprocessor, then we need two output
537  // slots for each input variable.
538  offset += this->dof_data[dataset]->n_output_variables *
539  (this->dof_data[dataset]->is_complex_valued() &&
540  (this->dof_data[dataset]->postprocessor == nullptr) ?
541  2 :
542  1);
543  }
544 
545  // then do the cell data. only compute the number of a cell if needed;
546  // also make sure that we only access cell data if the
547  // first_cell/next_cell functions only return active cells
548  if (this->cell_data.size() != 0)
549  {
550  Assert(!cell_and_index->first->has_children(), ExcNotImplemented());
551 
552  for (unsigned int dataset = 0; dataset < this->cell_data.size();
553  ++dataset)
554  {
555  // as above, first output the real part
556  {
557  const double value =
558  this->cell_data[dataset]->get_cell_data_value(
559  cell_and_index->second,
560  internal::DataOutImplementation::ComponentExtractor::
561  real_part);
562  for (unsigned int q = 0; q < n_q_points; ++q)
563  patch.data(offset, q) = value;
564  }
565 
566  // and if there is one, also output the imaginary part
567  if (this->cell_data[dataset]->is_complex_valued() == true)
568  {
569  const double value =
570  this->cell_data[dataset]->get_cell_data_value(
571  cell_and_index->second,
572  internal::DataOutImplementation::ComponentExtractor::
573  imaginary_part);
574  for (unsigned int q = 0; q < n_q_points; ++q)
575  patch.data(offset + 1, q) = value;
576  }
577 
578  offset += (this->cell_data[dataset]->is_complex_valued() ? 2 : 1);
579  }
580  }
581  }
582 
583 
584  for (unsigned int f = 0;
585  f < GeometryInfo<DoFHandlerType::dimension>::faces_per_cell;
586  ++f)
587  {
588  // let's look up whether the neighbor behind that face is noted in the
589  // table of cells which we treat. this can only happen if the neighbor
590  // exists, and is on the same level as this cell, but it may also happen
591  // that the neighbor is not a member of the range of cells over which we
592  // loop, in which case the respective entry in the
593  // cell_to_patch_index_map will have the value no_neighbor. (note that
594  // since we allocated only as much space in this array as the maximum
595  // index of the cells we loop over, not every neighbor may have its
596  // space in it, so we have to assume that it is extended by values
597  // no_neighbor)
598  if (cell_and_index->first->at_boundary(f) ||
599  (cell_and_index->first->neighbor(f)->level() !=
600  cell_and_index->first->level()))
601  {
602  patch.neighbors[f] = numbers::invalid_unsigned_int;
603  continue;
604  }
605 
606  const cell_iterator neighbor = cell_and_index->first->neighbor(f);
607  Assert(static_cast<unsigned int>(neighbor->level()) <
608  scratch_data.cell_to_patch_index_map->size(),
609  ExcInternalError());
610  if ((static_cast<unsigned int>(neighbor->index()) >=
611  (*scratch_data.cell_to_patch_index_map)[neighbor->level()].size()) ||
612  ((*scratch_data.cell_to_patch_index_map)[neighbor->level()]
613  [neighbor->index()] ==
615  {
616  patch.neighbors[f] = numbers::invalid_unsigned_int;
617  continue;
618  }
619 
620  // now, there is a neighbor, so get its patch number and set it for the
621  // neighbor index
622  patch.neighbors[f] =
623  (*scratch_data
624  .cell_to_patch_index_map)[neighbor->level()][neighbor->index()];
625  }
626 
627  const unsigned int patch_idx =
628  (*scratch_data.cell_to_patch_index_map)[cell_and_index->first->level()]
629  [cell_and_index->first->index()];
630  // did we mess up the indices?
631  Assert(patch_idx < this->patches.size(), ExcInternalError());
632  patch.patch_index = patch_idx;
633 
634  // Put the patch into the patches vector. instead of copying the data,
635  // simply swap the contents to avoid the penalty of writing into another
636  // processor's memory
637  this->patches[patch_idx].swap(patch);
638 }
639 
640 
641 
642 template <int dim, typename DoFHandlerType>
643 void
644 DataOut<dim, DoFHandlerType>::build_patches(const unsigned int n_subdivisions)
645 {
646  build_patches(StaticMappingQ1<DoFHandlerType::dimension,
647  DoFHandlerType::space_dimension>::mapping,
648  n_subdivisions,
649  no_curved_cells);
650 }
651 
652 
653 
654 template <int dim, typename DoFHandlerType>
655 void
658  & mapping,
659  const unsigned int n_subdivisions_,
660  const CurvedCellRegion curved_region)
661 {
662  // Check consistency of redundant template parameter
663  Assert(dim == DoFHandlerType::dimension,
664  ExcDimensionMismatch(dim, DoFHandlerType::dimension));
665 
666  Assert(this->triangulation != nullptr,
668 
669  const unsigned int n_subdivisions =
670  (n_subdivisions_ != 0) ? n_subdivisions_ : this->default_subdivisions;
671  Assert(n_subdivisions >= 1,
673  n_subdivisions));
674 
675  this->validate_dataset_names();
676 
677  // First count the cells we want to create patches of. Also fill the object
678  // that maps the cell indices to the patch numbers, as this will be needed
679  // for generation of neighborship information.
680  // Note, there is a confusing mess of different indices here at play:
681  // - patch_index: the index of a patch in all_cells
682  // - cell->index: only unique on each level, used in cell_to_patch_index_map
683  // - active_index: index for a cell when counting from begin_active() using
684  // ++cell (identical to cell->active_cell_index())
685  // - cell_index: unique index of a cell counted using
686  // next_locally_owned_cell() starting from first_locally_owned_cell()
687  //
688  // It turns out that we create one patch for each selected cell, so
689  // patch_index==cell_index.
690  //
691  // Now construct the map such that
692  // cell_to_patch_index_map[cell->level][cell->index] = patch_index
693  std::vector<std::vector<unsigned int>> cell_to_patch_index_map;
694  cell_to_patch_index_map.resize(this->triangulation->n_levels());
695  for (unsigned int l = 0; l < this->triangulation->n_levels(); ++l)
696  {
697  // max_index is the largest cell->index on level l
698  unsigned int max_index = 0;
699  for (cell_iterator cell = first_locally_owned_cell();
700  cell != this->triangulation->end();
701  cell = next_locally_owned_cell(cell))
702  if (static_cast<unsigned int>(cell->level()) == l)
703  max_index =
704  std::max(max_index, static_cast<unsigned int>(cell->index()));
705 
706  cell_to_patch_index_map[l].resize(
707  max_index + 1,
709  DoFHandlerType::dimension,
710  DoFHandlerType::space_dimension>::no_neighbor);
711  }
712 
713  // will be all_cells[patch_index] = pair(cell, active_index)
714  std::vector<std::pair<cell_iterator, unsigned int>> all_cells;
715  {
716  // important: we need to compute the active_index of the cell in the range
717  // 0..n_active_cells() because this is where we need to look up cell
718  // data from (cell data vectors do not have the length distance computed by
719  // first_locally_owned_cell/next_locally_owned_cell because this might skip
720  // some values (FilteredIterator).
721  active_cell_iterator active_cell = this->triangulation->begin_active();
722  unsigned int active_index = 0;
723  cell_iterator cell = first_locally_owned_cell();
724  for (; cell != this->triangulation->end();
725  cell = next_locally_owned_cell(cell))
726  {
727  // move forward until active_cell points at the cell (cell) we are
728  // looking at to compute the current active_index
729  while (active_cell != this->triangulation->end() && cell->active() &&
730  active_cell_iterator(cell) != active_cell)
731  {
732  ++active_cell;
733  ++active_index;
734  }
735 
736  Assert(static_cast<unsigned int>(cell->level()) <
737  cell_to_patch_index_map.size(),
738  ExcInternalError());
739  Assert(static_cast<unsigned int>(cell->index()) <
740  cell_to_patch_index_map[cell->level()].size(),
741  ExcInternalError());
742  Assert(active_index < this->triangulation->n_active_cells(),
743  ExcInternalError());
744  cell_to_patch_index_map[cell->level()][cell->index()] =
745  all_cells.size();
746 
747  all_cells.emplace_back(cell, active_index);
748  }
749  }
750 
751  this->patches.clear();
752  this->patches.resize(all_cells.size());
753 
754  // Now create a default object for the WorkStream object to work with. The
755  // first step is to count how many output data sets there will be. This is,
756  // in principle, just the number of components of each data set, but we
757  // need to allocate two entries per component if there are
758  // complex-valued input data (unless we use a postprocessor on this
759  // output -- all postprocessor outputs are real-valued)
760  unsigned int n_datasets = 0;
761  for (unsigned int i = 0; i < this->cell_data.size(); ++i)
762  n_datasets += (this->cell_data[i]->is_complex_valued() &&
763  (this->cell_data[i]->postprocessor == nullptr) ?
764  2 :
765  1);
766  for (unsigned int i = 0; i < this->dof_data.size(); ++i)
767  n_datasets += (this->dof_data[i]->n_output_variables *
768  (this->dof_data[i]->is_complex_valued() &&
769  (this->dof_data[i]->postprocessor == nullptr) ?
770  2 :
771  1));
772 
773  std::vector<unsigned int> n_postprocessor_outputs(this->dof_data.size());
774  for (unsigned int dataset = 0; dataset < this->dof_data.size(); ++dataset)
775  if (this->dof_data[dataset]->postprocessor)
776  n_postprocessor_outputs[dataset] =
777  this->dof_data[dataset]->n_output_variables;
778  else
779  n_postprocessor_outputs[dataset] = 0;
780 
781  const CurvedCellRegion curved_cell_region =
782  (n_subdivisions < 2 ? no_curved_cells : curved_region);
783 
784  UpdateFlags update_flags = update_values;
785  if (curved_cell_region != no_curved_cells)
786  update_flags |= update_quadrature_points;
787 
788  for (unsigned int i = 0; i < this->dof_data.size(); ++i)
789  if (this->dof_data[i]->postprocessor)
790  update_flags |=
791  this->dof_data[i]->postprocessor->get_needed_update_flags();
792  // perhaps update_normal_vectors is present, which would only be useful on
793  // faces, but we may not use it here.
794  Assert(
795  !(update_flags & update_normal_vectors),
796  ExcMessage(
797  "The update of normal vectors may not be requested for evaluation of "
798  "data on cells via DataPostprocessor."));
799 
800  internal::DataOutImplementation::ParallelData<DoFHandlerType::dimension,
801  DoFHandlerType::space_dimension>
802  thread_data(n_datasets,
803  n_subdivisions,
804  n_postprocessor_outputs,
805  mapping,
806  this->get_fes(),
807  update_flags,
808  cell_to_patch_index_map);
809 
810  // now build the patches in parallel
811  if (all_cells.size() > 0)
813  all_cells.data(),
814  all_cells.data() + all_cells.size(),
816  this,
817  std::placeholders::_1,
818  std::placeholders::_2,
819  /* no std::placeholders::_3, since this function doesn't
820  actually need a copy data object -- it just writes everything
821  right into the output array */
822  n_subdivisions,
823  curved_cell_region),
824  // no copy-local-to-global function needed here
825  std::function<void(const int)>(),
826  thread_data,
827  /* dummy CopyData object = */ 0,
828  // experimenting shows that we can make things run a bit
829  // faster if we increase the number of cells we work on
830  // per item (i.e., WorkStream's chunk_size argument,
831  // about 10% improvement) and the items in flight at any
832  // given time (another 5% on the testcase discussed in
833  // @ref workstream_paper, on 32 cores) and if
835  64);
836 }
837 
838 
839 
840 template <int dim, typename DoFHandlerType>
843 {
844  return this->triangulation->begin_active();
845 }
846 
847 
848 
849 template <int dim, typename DoFHandlerType>
852  const typename DataOut<dim, DoFHandlerType>::cell_iterator &cell)
853 {
854  // convert the iterator to an active_iterator and advance this to the next
855  // active cell
856  typename Triangulation<DoFHandlerType::dimension,
857  DoFHandlerType::space_dimension>::active_cell_iterator
858  active_cell = cell;
859  ++active_cell;
860  return active_cell;
861 }
862 
863 
864 
865 template <int dim, typename DoFHandlerType>
868 {
869  typename DataOut<dim, DoFHandlerType>::cell_iterator cell = first_cell();
870 
871  // skip cells if the current one has no children (is active) and is a ghost
872  // or artificial cell
873  while ((cell != this->triangulation->end()) &&
874  (cell->has_children() == false) && !cell->is_locally_owned())
875  cell = next_cell(cell);
876 
877  return cell;
878 }
879 
880 
881 
882 template <int dim, typename DoFHandlerType>
885  const typename DataOut<dim, DoFHandlerType>::cell_iterator &old_cell)
886 {
888  next_cell(old_cell);
889  while ((cell != this->triangulation->end()) &&
890  (cell->has_children() == false) && !cell->is_locally_owned())
891  cell = next_cell(cell);
892  return cell;
893 }
894 
895 
896 // explicit instantiations
897 #include "data_out.inst"
898 
899 DEAL_II_NAMESPACE_CLOSE
Shape function values.
static const unsigned int invalid_unsigned_int
Definition: types.h:187
static ::ExceptionBase & ExcNoTriangulationSelected()
void build_one_patch(const std::pair< cell_iterator, unsigned int > *cell_and_index, internal::DataOutImplementation::ParallelData< DoFHandlerType::dimension, DoFHandlerType::space_dimension > &scratch_data, const unsigned int n_subdivisions, const CurvedCellRegion curved_cell_region)
Definition: data_out.cc:70
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
virtual cell_iterator next_cell(const cell_iterator &cell)
Definition: data_out.cc:851
typename DataOut_DoFData< DoFHandler< dim >, DoFHandler< dim > ::dimension, DoFHandler< dim > ::space_dimension >::cell_iterator cell_iterator
Definition: data_out.h:172
virtual void evaluate_vector_field(const DataPostprocessorInputs::Vector< dim > &input_data, std::vector< Vector< double >> &computed_quantities) const
Transformed quadrature points.
static ::ExceptionBase & ExcInvalidNumberOfSubdivisions(int arg1)
virtual void build_patches(const unsigned int n_subdivisions=0)
Definition: data_out.cc:644
const std::vector< Point< spacedim > > & get_quadrature_points() const
static ::ExceptionBase & ExcMessage(std::string arg1)
virtual cell_iterator first_locally_owned_cell()
Definition: data_out.cc:867
#define Assert(cond, exc)
Definition: exceptions.h:1407
UpdateFlags
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Abstract base class for mapping classes.
Definition: mapping.h:301
virtual cell_iterator first_cell()
Definition: data_out.cc:842
unsigned int n_subdivisions
Second derivatives of shape functions.
const unsigned int n_quadrature_points
Definition: fe_values.h:2102
Shape function gradients.
virtual void evaluate_scalar_field(const DataPostprocessorInputs::Scalar< dim > &input_data, std::vector< Vector< double >> &computed_quantities) const
Normal vectors.
static ::ExceptionBase & ExcNotImplemented()
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:1167
static unsigned int n_threads()
virtual cell_iterator next_locally_owned_cell(const cell_iterator &cell)
Definition: data_out.cc:884
static ::ExceptionBase & ExcInternalError()
virtual UpdateFlags get_needed_update_flags() const =0