Reference documentation for deal.II version GIT 2896a7e638 2023-05-31 13:10:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
data_out_rotation.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 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 
18 
21 
22 #include <deal.II/fe/fe.h>
23 #include <deal.II/fe/fe_values.h>
24 #include <deal.II/fe/mapping.h>
25 
26 #include <deal.II/grid/tria.h>
28 
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<
90  &new_patches,
91  std::vector<
94  &patches)
95  {
96  for (unsigned int i = 0; i < new_patches.size(); ++i)
97  {
98  patches.push_back(new_patches[i]);
99  patches.back().patch_index = patches.size() - 1;
100  }
101  }
102  } // namespace DataOutRotationImplementation
103 } // namespace internal
104 
105 
106 
107 template <int dim, int spacedim>
108 void
110  const cell_iterator * cell,
112  std::vector<DataOutBase::Patch<patch_dim, patch_spacedim>> &my_patches)
113 {
114  if (dim == 3)
115  {
116  // would this function make any sense after all? who would want to
117  // output/compute in four space dimensions?
118  Assert(false, ExcNotImplemented());
119  return;
120  }
121 
122  Assert((*cell)->is_locally_owned(), ExcNotImplemented());
123 
124  const unsigned int n_patches_per_circle = data.n_patches_per_circle;
125 
126  // another abbreviation denoting the number of q_points in each direction
127  const unsigned int n_points = data.n_subdivisions + 1;
128 
129  // set up an array that holds the directions in the plane of rotation in
130  // which we will put points in the whole domain (not the rotationally
131  // reduced one in which the computation took place. for simplicity add the
132  // initial direction at the end again
133  std::vector<Point<dim + 1>> angle_directions(n_patches_per_circle + 1);
134  for (unsigned int i = 0; i <= n_patches_per_circle; ++i)
135  {
136  angle_directions[i][dim - 1] =
137  std::cos(2 * numbers::PI * i / n_patches_per_circle);
138  angle_directions[i][dim] =
139  std::sin(2 * numbers::PI * i / n_patches_per_circle);
140  }
141 
142  for (unsigned int angle = 0; angle < n_patches_per_circle; ++angle)
143  {
144  // first compute the vertices of the patch. note that they will have to
145  // be computed from the vertices of the cell, which has one dim
146  // less, however.
147  switch (dim)
148  {
149  case 1:
150  {
151  const double r1 = (*cell)->vertex(0)(0),
152  r2 = (*cell)->vertex(1)(0);
153  Assert(r1 >= 0, ExcRadialVariableHasNegativeValues(r1));
154  Assert(r2 >= 0, ExcRadialVariableHasNegativeValues(r2));
155 
156  my_patches[angle].vertices[0] = r1 * angle_directions[angle];
157  my_patches[angle].vertices[1] = r2 * angle_directions[angle];
158  my_patches[angle].vertices[2] = r1 * angle_directions[angle + 1];
159  my_patches[angle].vertices[3] = r2 * angle_directions[angle + 1];
160 
161  break;
162  }
163 
164  case 2:
165  {
166  for (const unsigned int vertex :
168  {
169  const Point<dim> v = (*cell)->vertex(vertex);
170 
171  // make sure that the radial variable is nonnegative
172  Assert(v(0) >= 0, ExcRadialVariableHasNegativeValues(v(0)));
173 
174  // now set the vertices of the patch
175  my_patches[angle].vertices[vertex] =
176  v(0) * angle_directions[angle];
177  my_patches[angle].vertices[vertex][0] = v(1);
178 
179  my_patches[angle]
180  .vertices[vertex + GeometryInfo<dim>::vertices_per_cell] =
181  v(0) * angle_directions[angle + 1];
182  my_patches[angle]
183  .vertices[vertex + GeometryInfo<dim>::vertices_per_cell]
184  [0] = v(1);
185  }
186 
187  break;
188  }
189 
190  default:
191  Assert(false, ExcNotImplemented());
192  }
193 
194  // then fill in data
195  if (data.n_datasets > 0)
196  {
197  unsigned int offset = 0;
198 
199  data.reinit_all_fe_values(this->dof_data, *cell);
200  // first fill dof_data
201  for (unsigned int dataset = 0; dataset < this->dof_data.size();
202  ++dataset)
203  {
204  const FEValuesBase<dim> &fe_patch_values =
205  data.get_present_fe_values(dataset);
206  const unsigned int n_components =
207  fe_patch_values.get_fe().n_components();
208  const DataPostprocessor<dim> *postprocessor =
209  this->dof_data[dataset]->postprocessor;
210  if (postprocessor != nullptr)
211  {
212  // we have to postprocess the
213  // data, so determine, which
214  // fields have to be updated
215  const UpdateFlags update_flags =
216  postprocessor->get_needed_update_flags();
217 
218  if (n_components == 1)
219  {
220  // at each point there is
221  // only one component of
222  // value, gradient etc.
223  if (update_flags & update_values)
224  this->dof_data[dataset]->get_function_values(
225  fe_patch_values,
227  real_part,
228  data.patch_values_scalar.solution_values);
229  if (update_flags & update_gradients)
230  this->dof_data[dataset]->get_function_gradients(
231  fe_patch_values,
233  real_part,
234  data.patch_values_scalar.solution_gradients);
235  if (update_flags & update_hessians)
236  this->dof_data[dataset]->get_function_hessians(
237  fe_patch_values,
239  real_part,
240  data.patch_values_scalar.solution_hessians);
241 
242  if (update_flags & update_quadrature_points)
243  data.patch_values_scalar.evaluation_points =
244  fe_patch_values.get_quadrature_points();
245 
246  const typename DoFHandler<dim,
247  spacedim>::active_cell_iterator
248  dh_cell(&(*cell)->get_triangulation(),
249  (*cell)->level(),
250  (*cell)->index(),
251  this->dof_data[dataset]->dof_handler);
252  data.patch_values_scalar.template set_cell<dim>(dh_cell);
253 
254  postprocessor->evaluate_scalar_field(
255  data.patch_values_scalar,
256  data.postprocessed_values[dataset]);
257  }
258  else
259  {
260  data.resize_system_vectors(n_components);
261 
262  // at each point there is a vector valued function and
263  // its derivative...
264  if (update_flags & update_values)
265  this->dof_data[dataset]->get_function_values(
266  fe_patch_values,
268  real_part,
269  data.patch_values_system.solution_values);
270  if (update_flags & update_gradients)
271  this->dof_data[dataset]->get_function_gradients(
272  fe_patch_values,
274  real_part,
275  data.patch_values_system.solution_gradients);
276  if (update_flags & update_hessians)
277  this->dof_data[dataset]->get_function_hessians(
278  fe_patch_values,
280  real_part,
281  data.patch_values_system.solution_hessians);
282 
283  if (update_flags & update_quadrature_points)
284  data.patch_values_system.evaluation_points =
285  fe_patch_values.get_quadrature_points();
286 
287  const typename DoFHandler<dim,
288  spacedim>::active_cell_iterator
289  dh_cell(&(*cell)->get_triangulation(),
290  (*cell)->level(),
291  (*cell)->index(),
292  this->dof_data[dataset]->dof_handler);
293  data.patch_values_system.template set_cell<dim>(dh_cell);
294 
295  postprocessor->evaluate_vector_field(
296  data.patch_values_system,
297  data.postprocessed_values[dataset]);
298  }
299 
300  for (unsigned int component = 0;
301  component < this->dof_data[dataset]->n_output_variables;
302  ++component)
303  {
304  switch (dim)
305  {
306  case 1:
307  for (unsigned int x = 0; x < n_points; ++x)
308  for (unsigned int y = 0; y < n_points; ++y)
309  my_patches[angle].data(offset + component,
310  x * n_points + y) =
311  data.postprocessed_values[dataset][x](
312  component);
313  break;
314 
315  case 2:
316  for (unsigned int x = 0; x < n_points; ++x)
317  for (unsigned int y = 0; y < n_points; ++y)
318  for (unsigned int z = 0; z < n_points; ++z)
319  my_patches[angle].data(offset + component,
320  x * n_points *
321  n_points +
322  y * n_points + z) =
323  data.postprocessed_values[dataset]
324  [x * n_points + z](
325  component);
326  break;
327 
328  default:
329  Assert(false, ExcNotImplemented());
330  }
331  }
332  }
333  else if (n_components == 1)
334  {
335  this->dof_data[dataset]->get_function_values(
336  fe_patch_values,
338  real_part,
339  data.patch_values_scalar.solution_values);
340 
341  switch (dim)
342  {
343  case 1:
344  for (unsigned int x = 0; x < n_points; ++x)
345  for (unsigned int y = 0; y < n_points; ++y)
346  my_patches[angle].data(offset, x * n_points + y) =
347  data.patch_values_scalar.solution_values[x];
348  break;
349 
350  case 2:
351  for (unsigned int x = 0; x < n_points; ++x)
352  for (unsigned int y = 0; y < n_points; ++y)
353  for (unsigned int z = 0; z < n_points; ++z)
354  my_patches[angle].data(offset,
355  x * n_points * n_points +
356  y + z * n_points) =
358  .solution_values[x * n_points + z];
359  break;
360 
361  default:
362  Assert(false, ExcNotImplemented());
363  }
364  }
365  else
366  // system of components
367  {
368  data.resize_system_vectors(n_components);
369  this->dof_data[dataset]->get_function_values(
370  fe_patch_values,
372  real_part,
373  data.patch_values_system.solution_values);
374 
375  for (unsigned int component = 0; component < n_components;
376  ++component)
377  {
378  switch (dim)
379  {
380  case 1:
381  for (unsigned int x = 0; x < n_points; ++x)
382  for (unsigned int y = 0; y < n_points; ++y)
383  my_patches[angle].data(offset + component,
384  x * n_points + y) =
385  data.patch_values_system.solution_values[x](
386  component);
387  break;
388 
389  case 2:
390  for (unsigned int x = 0; x < n_points; ++x)
391  for (unsigned int y = 0; y < n_points; ++y)
392  for (unsigned int z = 0; z < n_points; ++z)
393  my_patches[angle].data(offset + component,
394  x * n_points *
395  n_points +
396  y * n_points + z) =
398  .solution_values[x * n_points + z](
399  component);
400  break;
401 
402  default:
403  Assert(false, ExcNotImplemented());
404  }
405  }
406  }
407  offset += this->dof_data[dataset]->n_output_variables;
408  }
409 
410  // then do the cell data
411  for (unsigned int dataset = 0; dataset < this->cell_data.size();
412  ++dataset)
413  {
414  // we need to get at the number of the cell to which this face
415  // belongs in order to access the cell data. this is not readily
416  // available, so choose the following rather inefficient way:
417  Assert((*cell)->is_active(),
418  ExcMessage("Cell must be active for cell data"));
419  const unsigned int cell_number = std::distance(
420  this->triangulation->begin_active(),
422  *cell));
423  const double value =
424  this->cell_data[dataset]->get_cell_data_value(
425  cell_number,
427  real_part);
428  switch (dim)
429  {
430  case 1:
431  for (unsigned int x = 0; x < n_points; ++x)
432  for (unsigned int y = 0; y < n_points; ++y)
433  my_patches[angle].data(dataset + offset,
434  x * n_points + y) = value;
435  break;
436 
437  case 2:
438  for (unsigned int x = 0; x < n_points; ++x)
439  for (unsigned int y = 0; y < n_points; ++y)
440  for (unsigned int z = 0; z < n_points; ++z)
441  my_patches[angle].data(dataset + offset,
442  x * n_points * n_points +
443  y * n_points + z) = value;
444  break;
445 
446  default:
447  Assert(false, ExcNotImplemented());
448  }
449  }
450  }
451  }
452 }
453 
454 
455 
456 template <int dim, int spacedim>
457 void
459  const unsigned int n_patches_per_circle,
460  const unsigned int nnnn_subdivisions)
461 {
462  Assert(this->triangulation != nullptr,
464 
465  const unsigned int n_subdivisions =
466  (nnnn_subdivisions != 0) ? nnnn_subdivisions : this->default_subdivisions;
467  Assert(n_subdivisions >= 1,
469  n_subdivisions));
470 
471  this->validate_dataset_names();
472 
473  unsigned int n_datasets = this->cell_data.size();
474  for (unsigned int i = 0; i < this->dof_data.size(); ++i)
475  n_datasets += this->dof_data[i]->n_output_variables;
476 
478  for (unsigned int i = 0; i < this->dof_data.size(); ++i)
479  if (this->dof_data[i]->postprocessor)
480  update_flags |=
481  this->dof_data[i]->postprocessor->get_needed_update_flags();
482  // perhaps update_normal_vectors is present,
483  // which would only be useful on faces, but
484  // we may not use it here.
485  Assert(!(update_flags & update_normal_vectors),
486  ExcMessage("The update of normal vectors may not be requested for "
487  "evaluation of data on cells via DataPostprocessor."));
488 
489  // first count the cells we want to
490  // create patches of and make sure
491  // there is enough memory for that
492  std::vector<cell_iterator> all_cells;
493  for (cell_iterator cell = first_cell(); cell != this->triangulation->end();
494  cell = next_cell(cell))
495  all_cells.push_back(cell);
496 
497  // then also take into account that
498  // we want more than one patch to
499  // come out of every cell, as they
500  // are repeated around the axis of
501  // rotation
502  this->patches.clear();
503  this->patches.reserve(all_cells.size() * n_patches_per_circle);
504 
505 
506  std::vector<unsigned int> n_postprocessor_outputs(this->dof_data.size());
507  for (unsigned int dataset = 0; dataset < this->dof_data.size(); ++dataset)
508  if (this->dof_data[dataset]->postprocessor)
509  n_postprocessor_outputs[dataset] =
510  this->dof_data[dataset]->n_output_variables;
511  else
512  n_postprocessor_outputs[dataset] = 0;
513 
514  Assert(!this->triangulation->is_mixed_mesh(), ExcNotImplemented());
515  const auto reference_cell = this->triangulation->get_reference_cells()[0];
517  thread_data(
518  n_datasets,
519  n_subdivisions,
520  n_patches_per_circle,
521  n_postprocessor_outputs,
522  reference_cell.template get_default_linear_mapping<dim, spacedim>(),
523  this->get_fes(),
524  update_flags);
525  std::vector<DataOutBase::Patch<patch_dim, patch_spacedim>> new_patches(
526  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].reference_cell = ReferenceCells::get_hypercube<dim + 1>();
531 
532  new_patches[i].data.reinit(
533  n_datasets, Utilities::fixed_power<patch_dim>(n_subdivisions + 1));
534  }
535 
536  // now build the patches in parallel
538  all_cells.data(),
539  all_cells.data() + all_cells.size(),
540  [this](
541  const cell_iterator *cell,
543  & data,
544  std::vector<DataOutBase::Patch<patch_dim, patch_spacedim>> &my_patches) {
545  this->build_one_patch(cell, data, my_patches);
546  },
547  [this](const std::vector<DataOutBase::Patch<patch_dim, patch_spacedim>>
548  &new_patches) {
549  internal::DataOutRotationImplementation::append_patch_to_list<dim,
550  spacedim>(
551  new_patches, this->patches);
552  },
553  thread_data,
554  new_patches);
555 }
556 
557 
558 
559 template <int dim, int spacedim>
562 {
563  return this->triangulation->begin_active();
564 }
565 
566 
567 template <int dim, int spacedim>
570 {
571  // convert the iterator to an
572  // active_iterator and advance
573  // this to the next active cell
575  cell;
576  ++active_cell;
577  return active_cell;
578 }
579 
580 
581 
582 // explicit instantiations
583 #include "data_out_rotation.inst"
584 
585 
virtual void build_patches(const unsigned int n_patches_per_circle, const unsigned int n_subdivisions=0)
typename DataOut_DoFData< dim, patch_dim, spacedim, patch_spacedim >::cell_iterator cell_iterator
virtual cell_iterator next_cell(const cell_iterator &cell)
virtual cell_iterator first_cell()
void build_one_patch(const cell_iterator *cell, internal::DataOutRotationImplementation::ParallelData< dim, spacedim > &data, std::vector< DataOutBase::Patch< patch_dim, patch_spacedim >> &my_patches)
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 FiniteElement< dim, spacedim > & get_fe() const
const std::vector< Point< spacedim > > & get_quadrature_points() const
unsigned int n_components() const
Definition: point.h:112
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcInvalidNumberOfSubdivisions(int arg1)
static ::ExceptionBase & ExcNoTriangulationSelected()
#define Assert(cond, exc)
Definition: exceptions.h:1614
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcMessage(std::string arg1)
UpdateFlags
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
Number angle(const Tensor< 1, spacedim, Number > &a, const Tensor< 1, spacedim, Number > &b)
void run(const std::vector< std::vector< Iterator >> &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: work_stream.h:1270
void append_patch_to_list(const std::vector< DataOutBase::Patch< DataOutRotation< dim, spacedim >::patch_dim, DataOutRotation< dim, spacedim >::patch_spacedim >> &new_patches, std::vector< DataOutBase::Patch< DataOutRotation< dim, spacedim >::patch_dim, DataOutRotation< dim, spacedim >::patch_spacedim >> &patches)
static constexpr double PI
Definition: numbers.h:259
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)
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)