Reference documentation for deal.II version 9.4.1
\(\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\}}\)
Loading...
Searching...
No Matches
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>
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
52namespace 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
107template <int dim, int spacedim>
108void
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(
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(
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
456template <int dim, int spacedim>
457void
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
515 thread_data(n_datasets,
516 n_subdivisions,
517 n_patches_per_circle,
518 n_postprocessor_outputs,
520 this->get_fes(),
521 update_flags);
522 std::vector<DataOutBase::Patch<patch_dim, patch_spacedim>> new_patches(
523 n_patches_per_circle);
524 for (unsigned int i = 0; i < new_patches.size(); ++i)
525 {
526 new_patches[i].n_subdivisions = n_subdivisions;
527 new_patches[i].reference_cell = ReferenceCells::get_hypercube<dim + 1>();
528
529 new_patches[i].data.reinit(
530 n_datasets, Utilities::fixed_power<patch_dim>(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](
538 const cell_iterator *cell,
540 & data,
541 std::vector<DataOutBase::Patch<patch_dim, patch_spacedim>> &my_patches) {
542 this->build_one_patch(cell, data, my_patches);
543 },
544 [this](const std::vector<DataOutBase::Patch<patch_dim, patch_spacedim>>
545 &new_patches) {
546 internal::DataOutRotationImplementation::append_patch_to_list<dim,
547 spacedim>(
548 new_patches, this->patches);
549 },
550 thread_data,
551 new_patches);
552}
553
554
555
556template <int dim, int spacedim>
559{
560 return this->triangulation->begin_active();
561}
562
563
564template <int dim, int spacedim>
567{
568 // convert the iterator to an
569 // active_iterator and advance
570 // this to the next active cell
572 cell;
573 ++active_cell;
574 return active_cell;
575}
576
577
578
579// explicit instantiations
580#include "data_out_rotation.inst"
581
582
virtual void build_patches(const unsigned int n_patches_per_circle, const unsigned int n_subdivisions=0)
void build_one_patch(const cell_iterator *cell, internal::DataOutRotationImplementation::ParallelData< dim, spacedim > &data, std::vector< DataOutBase::Patch< patch_dim, patch_spacedim > > &my_patches)
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()
virtual UpdateFlags get_needed_update_flags() const =0
virtual void evaluate_vector_field(const DataPostprocessorInputs::Vector< dim > &input_data, std::vector< Vector< double > > &computed_quantities) const
virtual void evaluate_scalar_field(const DataPostprocessorInputs::Scalar< dim > &input_data, std::vector< Vector< double > > &computed_quantities) const
const std::vector< Point< spacedim > > & get_quadrature_points() const
const FiniteElement< dim, spacedim > & get_fe() const
unsigned int n_components() const
Abstract base class for mapping classes.
Definition: mapping.h:311
Definition: point.h:111
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
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.
static ::ExceptionBase & ExcInvalidNumberOfSubdivisions(int arg1)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNoTriangulationSelected()
static ::ExceptionBase & ExcMessage(std::string arg1)
void run(const std::vector< std::vector< Iterator > > &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: work_stream.h:1275
void 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:233
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
const FEValuesBase< dim, spacedim > & get_present_fe_values(const unsigned int dataset) const
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)
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)