Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
data_out_rotation.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 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/quadrature_lib.h>
17 #include <deal.II/base/work_stream.h>
18 
19 #include <deal.II/dofs/dof_accessor.h>
20 #include <deal.II/dofs/dof_handler.h>
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>
27 #include <deal.II/grid/tria_iterator.h>
28 
29 #include <deal.II/hp/fe_values.h>
30 
31 #include <deal.II/lac/block_vector.h>
32 #include <deal.II/lac/vector.h>
33 
34 #include <deal.II/numerics/data_out_rotation.h>
35 
36 #include <cmath>
37 
38 DEAL_II_NAMESPACE_OPEN
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>
57  ParallelData<dim, spacedim>::ParallelData(
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
86  append_patch_to_list(
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 (unsigned int vertex = 0;
163  vertex < GeometryInfo<dimension>::vertices_per_cell;
164  ++vertex)
165  {
166  const Point<dimension> v = (*cell)->vertex(vertex);
167 
168  // make sure that the radial variable is nonnegative
169  Assert(v(0) >= 0, ExcRadialVariableHasNegativeValues(v(0)));
170 
171  // now set the vertices of the patch
172  my_patches[angle].vertices[vertex] =
173  v(0) * angle_directions[angle];
174  my_patches[angle].vertices[vertex][0] = v(1);
175 
176  my_patches[angle]
177  .vertices[vertex +
179  v(0) * angle_directions[angle + 1];
180  my_patches[angle]
181  .vertices[vertex +
183  v(1);
184  }
185 
186  break;
187  }
188 
189  default:
190  Assert(false, ExcNotImplemented());
191  }
192 
193  // then fill in data
194  if (data.n_datasets > 0)
195  {
196  unsigned int offset = 0;
197 
198  data.reinit_all_fe_values(this->dof_data, *cell);
199  // first fill dof_data
200  for (unsigned int dataset = 0; dataset < this->dof_data.size();
201  ++dataset)
202  {
203  const FEValuesBase<dimension> &fe_patch_values =
204  data.get_present_fe_values(dataset);
205  const unsigned int n_components =
206  fe_patch_values.get_fe().n_components();
207  const DataPostprocessor<dim> *postprocessor =
208  this->dof_data[dataset]->postprocessor;
209  if (postprocessor != nullptr)
210  {
211  // we have to postprocess the
212  // data, so determine, which
213  // fields have to be updated
214  const UpdateFlags update_flags =
215  postprocessor->get_needed_update_flags();
216 
217  if (n_components == 1)
218  {
219  // at each point there is
220  // only one component of
221  // value, gradient etc.
222  if (update_flags & update_values)
223  this->dof_data[dataset]->get_function_values(
224  fe_patch_values,
225  internal::DataOutImplementation::ComponentExtractor::
226  real_part,
227  data.patch_values_scalar.solution_values);
228  if (update_flags & update_gradients)
229  this->dof_data[dataset]->get_function_gradients(
230  fe_patch_values,
231  internal::DataOutImplementation::ComponentExtractor::
232  real_part,
233  data.patch_values_scalar.solution_gradients);
234  if (update_flags & update_hessians)
235  this->dof_data[dataset]->get_function_hessians(
236  fe_patch_values,
237  internal::DataOutImplementation::ComponentExtractor::
238  real_part,
239  data.patch_values_scalar.solution_hessians);
240 
241  if (update_flags & update_quadrature_points)
242  data.patch_values_scalar.evaluation_points =
243  fe_patch_values.get_quadrature_points();
244 
245  const typename DoFHandlerType::active_cell_iterator
246  dh_cell(&(*cell)->get_triangulation(),
247  (*cell)->level(),
248  (*cell)->index(),
249  this->dof_data[dataset]->dof_handler);
250  data.patch_values_scalar
251  .template set_cell<DoFHandlerType>(dh_cell);
252 
253  postprocessor->evaluate_scalar_field(
254  data.patch_values_scalar,
255  data.postprocessed_values[dataset]);
256  }
257  else
258  {
259  data.resize_system_vectors(n_components);
260 
261  // at each point there is a vector valued function and
262  // its derivative...
263  if (update_flags & update_values)
264  this->dof_data[dataset]->get_function_values(
265  fe_patch_values,
266  internal::DataOutImplementation::ComponentExtractor::
267  real_part,
268  data.patch_values_system.solution_values);
269  if (update_flags & update_gradients)
270  this->dof_data[dataset]->get_function_gradients(
271  fe_patch_values,
272  internal::DataOutImplementation::ComponentExtractor::
273  real_part,
274  data.patch_values_system.solution_gradients);
275  if (update_flags & update_hessians)
276  this->dof_data[dataset]->get_function_hessians(
277  fe_patch_values,
278  internal::DataOutImplementation::ComponentExtractor::
279  real_part,
280  data.patch_values_system.solution_hessians);
281 
282  if (update_flags & update_quadrature_points)
283  data.patch_values_system.evaluation_points =
284  fe_patch_values.get_quadrature_points();
285 
286  const typename DoFHandlerType::active_cell_iterator
287  dh_cell(&(*cell)->get_triangulation(),
288  (*cell)->level(),
289  (*cell)->index(),
290  this->dof_data[dataset]->dof_handler);
291  data.patch_values_system
292  .template set_cell<DoFHandlerType>(dh_cell);
293 
294  postprocessor->evaluate_vector_field(
295  data.patch_values_system,
296  data.postprocessed_values[dataset]);
297  }
298 
299  for (unsigned int component = 0;
300  component < this->dof_data[dataset]->n_output_variables;
301  ++component)
302  {
303  switch (dimension)
304  {
305  case 1:
306  for (unsigned int x = 0; x < n_points; ++x)
307  for (unsigned int y = 0; y < n_points; ++y)
308  my_patches[angle].data(offset + component,
309  x * n_points + y) =
310  data.postprocessed_values[dataset][x](
311  component);
312  break;
313 
314  case 2:
315  for (unsigned int x = 0; x < n_points; ++x)
316  for (unsigned int y = 0; y < n_points; ++y)
317  for (unsigned int z = 0; z < n_points; ++z)
318  my_patches[angle].data(offset + component,
319  x * n_points *
320  n_points +
321  y * n_points + z) =
322  data.postprocessed_values[dataset]
323  [x * n_points + z](
324  component);
325  break;
326 
327  default:
328  Assert(false, ExcNotImplemented());
329  }
330  }
331  }
332  else if (n_components == 1)
333  {
334  this->dof_data[dataset]->get_function_values(
335  fe_patch_values,
336  internal::DataOutImplementation::ComponentExtractor::
337  real_part,
338  data.patch_values_scalar.solution_values);
339 
340  switch (dimension)
341  {
342  case 1:
343  for (unsigned int x = 0; x < n_points; ++x)
344  for (unsigned int y = 0; y < n_points; ++y)
345  my_patches[angle].data(offset, x * n_points + y) =
346  data.patch_values_scalar.solution_values[x];
347  break;
348 
349  case 2:
350  for (unsigned int x = 0; x < n_points; ++x)
351  for (unsigned int y = 0; y < n_points; ++y)
352  for (unsigned int z = 0; z < n_points; ++z)
353  my_patches[angle].data(offset,
354  x * n_points * n_points +
355  y + z * n_points) =
356  data.patch_values_scalar
357  .solution_values[x * n_points + z];
358  break;
359 
360  default:
361  Assert(false, ExcNotImplemented());
362  }
363  }
364  else
365  // system of components
366  {
367  data.resize_system_vectors(n_components);
368  this->dof_data[dataset]->get_function_values(
369  fe_patch_values,
370  internal::DataOutImplementation::ComponentExtractor::
371  real_part,
372  data.patch_values_system.solution_values);
373 
374  for (unsigned int component = 0; component < n_components;
375  ++component)
376  {
377  switch (dimension)
378  {
379  case 1:
380  for (unsigned int x = 0; x < n_points; ++x)
381  for (unsigned int y = 0; y < n_points; ++y)
382  my_patches[angle].data(offset + component,
383  x * n_points + y) =
384  data.patch_values_system.solution_values[x](
385  component);
386  break;
387 
388  case 2:
389  for (unsigned int x = 0; x < n_points; ++x)
390  for (unsigned int y = 0; y < n_points; ++y)
391  for (unsigned int z = 0; z < n_points; ++z)
392  my_patches[angle].data(offset + component,
393  x * n_points *
394  n_points +
395  y * n_points + z) =
396  data.patch_values_system
397  .solution_values[x * n_points + z](
398  component);
399  break;
400 
401  default:
402  Assert(false, ExcNotImplemented());
403  }
404  }
405  }
406  offset += this->dof_data[dataset]->n_output_variables;
407  }
408 
409  // then do the cell data
410  for (unsigned int dataset = 0; dataset < this->cell_data.size();
411  ++dataset)
412  {
413  // we need to get at the number of the cell to which this face
414  // belongs in order to access the cell data. this is not readily
415  // available, so choose the following rather inefficient way:
416  Assert((*cell)->active(),
417  ExcMessage("Cell must be active for cell data"));
418  const unsigned int cell_number = std::distance(
419  this->triangulation->begin_active(),
421  active_cell_iterator(*cell));
422  const double value =
423  this->cell_data[dataset]->get_cell_data_value(
424  cell_number,
425  internal::DataOutImplementation::ComponentExtractor::
426  real_part);
427  switch (dimension)
428  {
429  case 1:
430  for (unsigned int x = 0; x < n_points; ++x)
431  for (unsigned int y = 0; y < n_points; ++y)
432  my_patches[angle].data(dataset + offset,
433  x * n_points + y) = value;
434  break;
435 
436  case 2:
437  for (unsigned int x = 0; x < n_points; ++x)
438  for (unsigned int y = 0; y < n_points; ++y)
439  for (unsigned int z = 0; z < n_points; ++z)
440  my_patches[angle].data(dataset + offset,
441  x * n_points * n_points +
442  y * n_points + z) = value;
443  break;
444 
445  default:
446  Assert(false, ExcNotImplemented());
447  }
448  }
449  }
450  }
451 }
452 
453 
454 
455 template <int dim, typename DoFHandlerType>
456 void
458  const unsigned int n_patches_per_circle,
459  const unsigned int nnnn_subdivisions)
460 {
461  // Check consistency of redundant
462  // template parameter
463  Assert(dim == dimension, ExcDimensionMismatch(dim, dimension));
464  Assert(this->triangulation != nullptr,
466 
467  const unsigned int n_subdivisions =
468  (nnnn_subdivisions != 0) ? nnnn_subdivisions : this->default_subdivisions;
469  Assert(n_subdivisions >= 1,
471  n_subdivisions));
472 
473  this->validate_dataset_names();
474 
475  unsigned int n_datasets = this->cell_data.size();
476  for (unsigned int i = 0; i < this->dof_data.size(); ++i)
477  n_datasets += this->dof_data[i]->n_output_variables;
478 
480  for (unsigned int i = 0; i < this->dof_data.size(); ++i)
481  if (this->dof_data[i]->postprocessor)
482  update_flags |=
483  this->dof_data[i]->postprocessor->get_needed_update_flags();
484  // perhaps update_normal_vectors is present,
485  // which would only be useful on faces, but
486  // we may not use it here.
487  Assert(!(update_flags & update_normal_vectors),
488  ExcMessage("The update of normal vectors may not be requested for "
489  "evaluation of data on cells via DataPostprocessor."));
490 
491  // first count the cells we want to
492  // create patches of and make sure
493  // there is enough memory for that
494  std::vector<cell_iterator> all_cells;
495  for (cell_iterator cell = first_cell(); cell != this->triangulation->end();
496  cell = next_cell(cell))
497  all_cells.push_back(cell);
498 
499  // then also take into account that
500  // we want more than one patch to
501  // come out of every cell, as they
502  // are repeated around the axis of
503  // rotation
504  this->patches.clear();
505  this->patches.reserve(all_cells.size() * n_patches_per_circle);
506 
507 
508  std::vector<unsigned int> n_postprocessor_outputs(this->dof_data.size());
509  for (unsigned int dataset = 0; dataset < this->dof_data.size(); ++dataset)
510  if (this->dof_data[dataset]->postprocessor)
511  n_postprocessor_outputs[dataset] =
512  this->dof_data[dataset]->n_output_variables;
513  else
514  n_postprocessor_outputs[dataset] = 0;
515 
517  space_dimension>
518  thread_data(n_datasets,
519  n_subdivisions,
520  n_patches_per_circle,
521  n_postprocessor_outputs,
523  this->get_fes(),
524  update_flags);
525  std::vector<DataOutBase::Patch<dimension + 1, space_dimension + 1>>
526  new_patches(n_patches_per_circle);
527  for (unsigned int i = 0; i < new_patches.size(); ++i)
528  {
529  new_patches[i].n_subdivisions = n_subdivisions;
530  new_patches[i].data.reinit(
531  n_datasets, Utilities::fixed_power<dimension + 1>(n_subdivisions + 1));
532  }
533 
534  // now build the patches in parallel
536  all_cells.data(),
537  all_cells.data() + all_cells.size(),
539  this,
540  std::placeholders::_1,
541  std::placeholders::_2,
542  std::placeholders::_3),
544  append_patch_to_list<dim, space_dimension>,
545  std::placeholders::_1,
546  std::ref(this->patches)),
547  thread_data,
548  new_patches);
549 }
550 
551 
552 
553 template <int dim, typename DoFHandlerType>
556 {
557  return this->triangulation->begin_active();
558 }
559 
560 
561 template <int dim, typename DoFHandlerType>
564 {
565  // convert the iterator to an
566  // active_iterator and advance
567  // this to the next active cell
569  active_cell = cell;
570  ++active_cell;
571  return active_cell;
572 }
573 
574 
575 
576 // explicit instantiations
577 #include "data_out_rotation.inst"
578 
579 
580 DEAL_II_NAMESPACE_CLOSE
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:1529
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:110
const std::vector< Point< spacedim > > & get_quadrature_points() const
static ::ExceptionBase & ExcMessage(std::string arg1)
#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: dof_tools.h:57
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:146
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:1167
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
virtual UpdateFlags get_needed_update_flags() const =0