Reference documentation for deal.II version 9.2.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_rotation.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 2020 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 
18 
21 
22 #include <deal.II/fe/fe.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/fe_values.h>
30 
32 #include <deal.II/lac/vector.h>
33 
35 
36 #include <cmath>
37 
39 
40 
41 // TODO: Update documentation
42 // TODO: Unify code for dimensions
43 
44 
45 // TODO: build_some_patches isn't going to work if first_cell/next_cell
46 // don't iterate over all cells and if cell data is requested. in that
47 // case, we need to calculate cell_number as in the DataOut class
48 
49 // Not implemented for 3D
50 
51 
52 namespace internal
53 {
54  namespace DataOutRotationImplementation
55  {
56  template <int dim, int spacedim>
58  const unsigned int n_datasets,
59  const unsigned int n_subdivisions,
60  const unsigned int n_patches_per_circle,
61  const std::vector<unsigned int> &n_postprocessor_outputs,
62  const Mapping<dim, spacedim> & mapping,
63  const std::vector<
64  std::shared_ptr<::hp::FECollection<dim, spacedim>>>
65  & finite_elements,
66  const UpdateFlags update_flags)
67  : internal::DataOutImplementation::ParallelDataBase<dim, spacedim>(
68  n_datasets,
69  n_subdivisions,
70  n_postprocessor_outputs,
71  mapping,
72  finite_elements,
73  update_flags,
74  false)
75  , n_patches_per_circle(n_patches_per_circle)
76  {}
77 
78 
79 
84  template <int dim, int spacedim>
85  void
87  const std::vector<DataOutBase::Patch<dim + 1, spacedim + 1>> &new_patches,
88  std::vector<DataOutBase::Patch<dim + 1, spacedim + 1>> & patches)
89  {
90  for (unsigned int i = 0; i < new_patches.size(); ++i)
91  {
92  patches.push_back(new_patches[i]);
93  patches.back().patch_index = patches.size() - 1;
94  }
95  }
96  } // namespace DataOutRotationImplementation
97 } // namespace internal
98 
99 
100 
101 template <int dim, typename DoFHandlerType>
102 void
104  const cell_iterator * cell,
106  space_dimension> &data,
108  &my_patches)
109 {
110  if (dim == 3)
111  {
112  // would this function make any sense after all? who would want to
113  // output/compute in four space dimensions?
114  Assert(false, ExcNotImplemented());
115  return;
116  }
117 
118  Assert((*cell)->is_locally_owned(), ExcNotImplemented());
119 
120  const unsigned int n_patches_per_circle = data.n_patches_per_circle;
121 
122  // another abbreviation denoting the number of q_points in each direction
123  const unsigned int n_points = data.n_subdivisions + 1;
124 
125  // set up an array that holds the directions in the plane of rotation in
126  // which we will put points in the whole domain (not the rotationally
127  // reduced one in which the computation took place. for simplicity add the
128  // initial direction at the end again
129  std::vector<Point<dimension + 1>> angle_directions(n_patches_per_circle + 1);
130  for (unsigned int i = 0; i <= n_patches_per_circle; ++i)
131  {
132  angle_directions[i][dimension - 1] =
133  std::cos(2 * numbers::PI * i / n_patches_per_circle);
134  angle_directions[i][dimension] =
135  std::sin(2 * numbers::PI * i / n_patches_per_circle);
136  }
137 
138  for (unsigned int angle = 0; angle < n_patches_per_circle; ++angle)
139  {
140  // first compute the vertices of the patch. note that they will have to
141  // be computed from the vertices of the cell, which has one dimension
142  // less, however.
143  switch (dimension)
144  {
145  case 1:
146  {
147  const double r1 = (*cell)->vertex(0)(0),
148  r2 = (*cell)->vertex(1)(0);
149  Assert(r1 >= 0, ExcRadialVariableHasNegativeValues(r1));
150  Assert(r2 >= 0, ExcRadialVariableHasNegativeValues(r2));
151 
152  my_patches[angle].vertices[0] = r1 * angle_directions[angle];
153  my_patches[angle].vertices[1] = r2 * angle_directions[angle];
154  my_patches[angle].vertices[2] = r1 * angle_directions[angle + 1];
155  my_patches[angle].vertices[3] = r2 * angle_directions[angle + 1];
156 
157  break;
158  }
159 
160  case 2:
161  {
162  for (const unsigned int vertex :
164  {
165  const Point<dimension> v = (*cell)->vertex(vertex);
166 
167  // make sure that the radial variable is nonnegative
168  Assert(v(0) >= 0, ExcRadialVariableHasNegativeValues(v(0)));
169 
170  // now set the vertices of the patch
171  my_patches[angle].vertices[vertex] =
172  v(0) * angle_directions[angle];
173  my_patches[angle].vertices[vertex][0] = v(1);
174 
175  my_patches[angle]
176  .vertices[vertex +
178  v(0) * angle_directions[angle + 1];
179  my_patches[angle]
180  .vertices[vertex +
182  v(1);
183  }
184 
185  break;
186  }
187 
188  default:
189  Assert(false, ExcNotImplemented());
190  }
191 
192  // then fill in data
193  if (data.n_datasets > 0)
194  {
195  unsigned int offset = 0;
196 
197  data.reinit_all_fe_values(this->dof_data, *cell);
198  // first fill dof_data
199  for (unsigned int dataset = 0; dataset < this->dof_data.size();
200  ++dataset)
201  {
202  const FEValuesBase<dimension> &fe_patch_values =
203  data.get_present_fe_values(dataset);
204  const unsigned int n_components =
205  fe_patch_values.get_fe().n_components();
206  const DataPostprocessor<dim> *postprocessor =
207  this->dof_data[dataset]->postprocessor;
208  if (postprocessor != nullptr)
209  {
210  // we have to postprocess the
211  // data, so determine, which
212  // fields have to be updated
213  const UpdateFlags update_flags =
214  postprocessor->get_needed_update_flags();
215 
216  if (n_components == 1)
217  {
218  // at each point there is
219  // only one component of
220  // value, gradient etc.
221  if (update_flags & update_values)
222  this->dof_data[dataset]->get_function_values(
223  fe_patch_values,
225  real_part,
226  data.patch_values_scalar.solution_values);
227  if (update_flags & update_gradients)
228  this->dof_data[dataset]->get_function_gradients(
229  fe_patch_values,
231  real_part,
232  data.patch_values_scalar.solution_gradients);
233  if (update_flags & update_hessians)
234  this->dof_data[dataset]->get_function_hessians(
235  fe_patch_values,
237  real_part,
238  data.patch_values_scalar.solution_hessians);
239 
240  if (update_flags & update_quadrature_points)
241  data.patch_values_scalar.evaluation_points =
242  fe_patch_values.get_quadrature_points();
243 
244  const typename DoFHandlerType::active_cell_iterator
245  dh_cell(&(*cell)->get_triangulation(),
246  (*cell)->level(),
247  (*cell)->index(),
248  this->dof_data[dataset]->dof_handler);
249  data.patch_values_scalar
250  .template set_cell<DoFHandlerType>(dh_cell);
251 
252  postprocessor->evaluate_scalar_field(
253  data.patch_values_scalar,
254  data.postprocessed_values[dataset]);
255  }
256  else
257  {
258  data.resize_system_vectors(n_components);
259 
260  // at each point there is a vector valued function and
261  // its derivative...
262  if (update_flags & update_values)
263  this->dof_data[dataset]->get_function_values(
264  fe_patch_values,
266  real_part,
267  data.patch_values_system.solution_values);
268  if (update_flags & update_gradients)
269  this->dof_data[dataset]->get_function_gradients(
270  fe_patch_values,
272  real_part,
273  data.patch_values_system.solution_gradients);
274  if (update_flags & update_hessians)
275  this->dof_data[dataset]->get_function_hessians(
276  fe_patch_values,
278  real_part,
279  data.patch_values_system.solution_hessians);
280 
281  if (update_flags & update_quadrature_points)
282  data.patch_values_system.evaluation_points =
283  fe_patch_values.get_quadrature_points();
284 
285  const typename DoFHandlerType::active_cell_iterator
286  dh_cell(&(*cell)->get_triangulation(),
287  (*cell)->level(),
288  (*cell)->index(),
289  this->dof_data[dataset]->dof_handler);
290  data.patch_values_system
291  .template set_cell<DoFHandlerType>(dh_cell);
292 
293  postprocessor->evaluate_vector_field(
294  data.patch_values_system,
295  data.postprocessed_values[dataset]);
296  }
297 
298  for (unsigned int component = 0;
299  component < this->dof_data[dataset]->n_output_variables;
300  ++component)
301  {
302  switch (dimension)
303  {
304  case 1:
305  for (unsigned int x = 0; x < n_points; ++x)
306  for (unsigned int y = 0; y < n_points; ++y)
307  my_patches[angle].data(offset + component,
308  x * n_points + y) =
309  data.postprocessed_values[dataset][x](
310  component);
311  break;
312 
313  case 2:
314  for (unsigned int x = 0; x < n_points; ++x)
315  for (unsigned int y = 0; y < n_points; ++y)
316  for (unsigned int z = 0; z < n_points; ++z)
317  my_patches[angle].data(offset + component,
318  x * n_points *
319  n_points +
320  y * n_points + z) =
321  data.postprocessed_values[dataset]
322  [x * n_points + z](
323  component);
324  break;
325 
326  default:
327  Assert(false, ExcNotImplemented());
328  }
329  }
330  }
331  else if (n_components == 1)
332  {
333  this->dof_data[dataset]->get_function_values(
334  fe_patch_values,
336  real_part,
337  data.patch_values_scalar.solution_values);
338 
339  switch (dimension)
340  {
341  case 1:
342  for (unsigned int x = 0; x < n_points; ++x)
343  for (unsigned int y = 0; y < n_points; ++y)
344  my_patches[angle].data(offset, x * n_points + y) =
345  data.patch_values_scalar.solution_values[x];
346  break;
347 
348  case 2:
349  for (unsigned int x = 0; x < n_points; ++x)
350  for (unsigned int y = 0; y < n_points; ++y)
351  for (unsigned int z = 0; z < n_points; ++z)
352  my_patches[angle].data(offset,
353  x * n_points * n_points +
354  y + z * n_points) =
355  data.patch_values_scalar
356  .solution_values[x * n_points + z];
357  break;
358 
359  default:
360  Assert(false, ExcNotImplemented());
361  }
362  }
363  else
364  // system of components
365  {
366  data.resize_system_vectors(n_components);
367  this->dof_data[dataset]->get_function_values(
368  fe_patch_values,
370  real_part,
371  data.patch_values_system.solution_values);
372 
373  for (unsigned int component = 0; component < n_components;
374  ++component)
375  {
376  switch (dimension)
377  {
378  case 1:
379  for (unsigned int x = 0; x < n_points; ++x)
380  for (unsigned int y = 0; y < n_points; ++y)
381  my_patches[angle].data(offset + component,
382  x * n_points + y) =
383  data.patch_values_system.solution_values[x](
384  component);
385  break;
386 
387  case 2:
388  for (unsigned int x = 0; x < n_points; ++x)
389  for (unsigned int y = 0; y < n_points; ++y)
390  for (unsigned int z = 0; z < n_points; ++z)
391  my_patches[angle].data(offset + component,
392  x * n_points *
393  n_points +
394  y * n_points + z) =
395  data.patch_values_system
396  .solution_values[x * n_points + z](
397  component);
398  break;
399 
400  default:
401  Assert(false, ExcNotImplemented());
402  }
403  }
404  }
405  offset += this->dof_data[dataset]->n_output_variables;
406  }
407 
408  // then do the cell data
409  for (unsigned int dataset = 0; dataset < this->cell_data.size();
410  ++dataset)
411  {
412  // we need to get at the number of the cell to which this face
413  // belongs in order to access the cell data. this is not readily
414  // available, so choose the following rather inefficient way:
415  Assert((*cell)->is_active(),
416  ExcMessage("Cell must be active for cell data"));
417  const unsigned int cell_number = std::distance(
418  this->triangulation->begin_active(),
420  active_cell_iterator(*cell));
421  const double value =
422  this->cell_data[dataset]->get_cell_data_value(
423  cell_number,
425  real_part);
426  switch (dimension)
427  {
428  case 1:
429  for (unsigned int x = 0; x < n_points; ++x)
430  for (unsigned int y = 0; y < n_points; ++y)
431  my_patches[angle].data(dataset + offset,
432  x * n_points + y) = value;
433  break;
434 
435  case 2:
436  for (unsigned int x = 0; x < n_points; ++x)
437  for (unsigned int y = 0; y < n_points; ++y)
438  for (unsigned int z = 0; z < n_points; ++z)
439  my_patches[angle].data(dataset + offset,
440  x * n_points * n_points +
441  y * n_points + z) = value;
442  break;
443 
444  default:
445  Assert(false, ExcNotImplemented());
446  }
447  }
448  }
449  }
450 }
451 
452 
453 
454 template <int dim, typename DoFHandlerType>
455 void
457  const unsigned int n_patches_per_circle,
458  const unsigned int nnnn_subdivisions)
459 {
460  // Check consistency of redundant
461  // template parameter
462  Assert(dim == dimension, ExcDimensionMismatch(dim, dimension));
463  Assert(this->triangulation != nullptr,
465 
466  const unsigned int n_subdivisions =
467  (nnnn_subdivisions != 0) ? nnnn_subdivisions : this->default_subdivisions;
468  Assert(n_subdivisions >= 1,
470  n_subdivisions));
471 
472  this->validate_dataset_names();
473 
474  unsigned int n_datasets = this->cell_data.size();
475  for (unsigned int i = 0; i < this->dof_data.size(); ++i)
476  n_datasets += this->dof_data[i]->n_output_variables;
477 
479  for (unsigned int i = 0; i < this->dof_data.size(); ++i)
480  if (this->dof_data[i]->postprocessor)
481  update_flags |=
482  this->dof_data[i]->postprocessor->get_needed_update_flags();
483  // perhaps update_normal_vectors is present,
484  // which would only be useful on faces, but
485  // we may not use it here.
486  Assert(!(update_flags & update_normal_vectors),
487  ExcMessage("The update of normal vectors may not be requested for "
488  "evaluation of data on cells via DataPostprocessor."));
489 
490  // first count the cells we want to
491  // create patches of and make sure
492  // there is enough memory for that
493  std::vector<cell_iterator> all_cells;
494  for (cell_iterator cell = first_cell(); cell != this->triangulation->end();
495  cell = next_cell(cell))
496  all_cells.push_back(cell);
497 
498  // then also take into account that
499  // we want more than one patch to
500  // come out of every cell, as they
501  // are repeated around the axis of
502  // rotation
503  this->patches.clear();
504  this->patches.reserve(all_cells.size() * n_patches_per_circle);
505 
506 
507  std::vector<unsigned int> n_postprocessor_outputs(this->dof_data.size());
508  for (unsigned int dataset = 0; dataset < this->dof_data.size(); ++dataset)
509  if (this->dof_data[dataset]->postprocessor)
510  n_postprocessor_outputs[dataset] =
511  this->dof_data[dataset]->n_output_variables;
512  else
513  n_postprocessor_outputs[dataset] = 0;
514 
516  space_dimension>
517  thread_data(n_datasets,
518  n_subdivisions,
519  n_patches_per_circle,
520  n_postprocessor_outputs,
522  this->get_fes(),
523  update_flags);
524  std::vector<DataOutBase::Patch<dimension + 1, space_dimension + 1>>
525  new_patches(n_patches_per_circle);
526  for (unsigned int i = 0; i < new_patches.size(); ++i)
527  {
528  new_patches[i].n_subdivisions = n_subdivisions;
529  new_patches[i].data.reinit(
530  n_datasets, Utilities::fixed_power<dimension + 1>(n_subdivisions + 1));
531  }
532 
533  // now build the patches in parallel
535  all_cells.data(),
536  all_cells.data() + all_cells.size(),
537  [this](const cell_iterator *cell,
538  internal::DataOutRotationImplementation::
539  ParallelData<dimension, space_dimension> &data,
540  std::vector<DataOutBase::Patch<dimension + 1, space_dimension + 1>>
541  &my_patches) { this->build_one_patch(cell, data, my_patches); },
542  [this](
543  const std::vector<DataOutBase::Patch<dimension + 1, space_dimension + 1>>
544  &new_patches) {
545  internal::DataOutRotationImplementation::
546  append_patch_to_list<dimension, space_dimension>(new_patches,
547  this->patches);
548  },
549  thread_data,
550  new_patches);
551 }
552 
553 
554 
555 template <int dim, typename DoFHandlerType>
558 {
559  return this->triangulation->begin_active();
560 }
561 
562 
563 template <int dim, typename DoFHandlerType>
566 {
567  // convert the iterator to an
568  // active_iterator and advance
569  // this to the next active cell
571  active_cell = cell;
572  ++active_cell;
573  return active_cell;
574 }
575 
576 
577 
578 // explicit instantiations
579 #include "data_out_rotation.inst"
580 
581 
Shape function values.
static ::ExceptionBase & ExcNoTriangulationSelected()
virtual cell_iterator next_cell(const cell_iterator &cell)
const FiniteElement< dim, spacedim > & get_fe() const
virtual void evaluate_vector_field(const DataPostprocessorInputs::Vector< dim > &input_data, std::vector< Vector< double >> &computed_quantities) const
Transformed quadrature points.
TriaActiveIterator< CellAccessor< dim, spacedim > > active_cell_iterator
Definition: tria.h:1361
void build_one_patch(const cell_iterator *cell, internal::DataOutRotationImplementation::ParallelData< dimension, space_dimension > &data, std::vector< DataOutBase::Patch< dimension+1, space_dimension+1 >> &my_patches)
static ::ExceptionBase & ExcInvalidNumberOfSubdivisions(int arg1)
Definition: point.h:111
const double angle
void append_patch_to_list(const std::vector< DataOutBase::Patch< dim+1, spacedim+1 >> &new_patches, std::vector< DataOutBase::Patch< dim+1, spacedim+1 >> &patches)
const std::vector< Point< spacedim > > & get_quadrature_points() const
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1419
UpdateFlags
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Abstract base class for mapping classes.
Definition: mapping.h:302
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
Second derivatives of shape functions.
virtual cell_iterator first_cell()
virtual void build_patches(const unsigned int n_patches_per_circle, const unsigned int n_subdivisions=0)
static constexpr double PI
Definition: numbers.h:237
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
Shape function gradients.
virtual void evaluate_scalar_field(const DataPostprocessorInputs::Scalar< dim > &input_data, std::vector< Vector< double >> &computed_quantities) const
Normal vectors.
typename DataOut_DoFData< DoFHandlerType, dimension+1 >::cell_iterator cell_iterator
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:1185
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
static const bool value
virtual UpdateFlags get_needed_update_flags() const =0
ParallelData(const unsigned int n_datasets, const unsigned int n_subdivisions, const unsigned int n_patches_per_circle, const std::vector< unsigned int > &n_postprocessor_outputs, const Mapping< dim, spacedim > &mapping, const std::vector< std::shared_ptr<::hp::FECollection< dim, spacedim >>> &finite_elements, const UpdateFlags update_flags)