deal.II version GIT relicensing-6842-g793a97d2aa 2026-10-02 14:00:01+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\}}\)
Loading...
Searching...
No Matches
data_out_rotation.cc
Go to the documentation of this file.
1// -----------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 2000 - 2026 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
15
18
19#include <deal.II/fe/fe.h>
21#include <deal.II/fe/mapping.h>
22
23#include <deal.II/grid/tria.h>
25
27
29#include <deal.II/lac/vector.h>
30
32
33#include <cmath>
34
36
37
38// TODO: Update documentation
39// TODO: Unify code for dimensions
40
41
42// TODO: build_some_patches isn't going to work if first_cell/next_cell
43// don't iterate over all cells and if cell data is requested. in that
44// case, we need to calculate cell_number as in the DataOut class
45
46// Not implemented for 3d
47
48
49namespace internal
50{
51 namespace DataOutRotationImplementation
52 {
53 template <int dim, int spacedim>
55 const unsigned int n_datasets,
56 const unsigned int n_subdivisions,
57 const unsigned int n_patches_per_circle,
58 const std::vector<unsigned int> &n_postprocessor_outputs,
59 const Mapping<dim, spacedim> &mapping,
60 const std::vector<
61 std::shared_ptr<::hp::FECollection<dim, spacedim>>>
62 &finite_elements,
63 const UpdateFlags update_flags)
64 : internal::DataOutImplementation::ParallelDataBase<dim, spacedim>(
65 n_datasets,
66 n_subdivisions,
67 n_postprocessor_outputs,
68 mapping,
69 finite_elements,
70 update_flags,
71 false)
72 , n_patches_per_circle(n_patches_per_circle)
73 {}
74
75
76
81 template <int dim, int spacedim>
82 void
84 const std::vector<
87 &new_patches,
88 std::vector<
91 &patches)
92 {
93 for (unsigned int i = 0; i < new_patches.size(); ++i)
94 {
95 patches.push_back(new_patches[i]);
96 patches.back().patch_index = patches.size() - 1;
97 }
98 }
99 } // namespace DataOutRotationImplementation
100} // namespace internal
101
102
103
104template <int dim, int spacedim>
105void
107 const cell_iterator *cell,
109 std::vector<DataOutBase::Patch<patch_dim, patch_spacedim>> &my_patches)
110{
111 if (dim == 3)
112 {
113 // would this function make any sense after all? who would want to
114 // output/compute in four space dimensions?
116 return;
117 }
118
119 Assert((*cell)->is_locally_owned(), ExcNotImplemented());
120
121 const unsigned int n_patches_per_circle = data.n_patches_per_circle;
122
123 // another abbreviation denoting the number of q_points in each direction
124 const unsigned int n_points = data.n_subdivisions + 1;
125
126 // set up an array that holds the directions in the plane of rotation in
127 // which we will put points in the whole domain (not the rotationally
128 // reduced one in which the computation took place. for simplicity add the
129 // initial direction at the end again
130 std::vector<Point<dim + 1>> angle_directions(n_patches_per_circle + 1);
131 for (unsigned int i = 0; i <= n_patches_per_circle; ++i)
132 {
133 angle_directions[i][dim - 1] =
134 std::cos(2 * numbers::PI * i / n_patches_per_circle);
135 angle_directions[i][dim] =
136 std::sin(2 * numbers::PI * i / n_patches_per_circle);
137 }
138
139 for (unsigned int angle = 0; angle < n_patches_per_circle; ++angle)
140 {
141 // first compute the vertices of the patch. note that they will have to
142 // be computed from the vertices of the cell, which has one dim
143 // less, however.
144 switch (dim)
145 {
146 case 1:
147 {
148 const double r1 = (*cell)->vertex(0)[0],
149 r2 = (*cell)->vertex(1)[0];
150 Assert(r1 >= 0, ExcRadialVariableHasNegativeValues(r1));
151 Assert(r2 >= 0, ExcRadialVariableHasNegativeValues(r2));
152
153 my_patches[angle].vertices[0] = r1 * angle_directions[angle];
154 my_patches[angle].vertices[1] = r2 * angle_directions[angle];
155 my_patches[angle].vertices[2] = r1 * angle_directions[angle + 1];
156 my_patches[angle].vertices[3] = r2 * angle_directions[angle + 1];
157
158 break;
159 }
160
161 case 2:
162 {
163 for (const unsigned int vertex :
165 {
166 const Point<dim> 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 + GeometryInfo<dim>::vertices_per_cell] =
178 v[0] * angle_directions[angle + 1];
179 my_patches[angle]
180 .vertices[vertex + GeometryInfo<dim>::vertices_per_cell]
181 [0] = v[1];
182 }
183
184 break;
185 }
186
187 default:
189 }
190
191 // then fill in data
192 if (data.n_datasets > 0)
193 {
194 unsigned int offset = 0;
195
196 data.reinit_all_fe_values(this->dof_data, *cell);
197 // first fill dof_data
198 for (unsigned int dataset = 0; dataset < this->dof_data.size();
199 ++dataset)
200 {
201 const FEValuesBase<dim> &fe_patch_values =
202 data.get_present_fe_values(dataset);
203 const unsigned int n_components =
204 fe_patch_values.get_fe().n_components();
205 const DataPostprocessor<dim> *postprocessor =
206 this->dof_data[dataset]->postprocessor;
207 if (postprocessor != nullptr)
208 {
209 // we have to postprocess the
210 // data, so determine, which
211 // fields have to be updated
212 const UpdateFlags update_flags =
213 postprocessor->get_needed_update_flags();
214
215 if (n_components == 1)
216 {
217 // at each point there is
218 // only one component of
219 // value, gradient etc.
220 if (update_flags & update_values)
221 this->dof_data[dataset]->get_function_values(
222 fe_patch_values,
224 real_part,
225 data.patch_values_scalar.solution_values);
226 if (update_flags & update_gradients)
227 this->dof_data[dataset]->get_function_gradients(
228 fe_patch_values,
230 real_part,
231 data.patch_values_scalar.solution_gradients);
232 if (update_flags & update_hessians)
233 this->dof_data[dataset]->get_function_hessians(
234 fe_patch_values,
236 real_part,
237 data.patch_values_scalar.solution_hessians);
238
239 if (update_flags & update_quadrature_points)
240 data.patch_values_scalar.evaluation_points =
241 fe_patch_values.get_quadrature_points();
242
243 const typename DoFHandler<dim,
244 spacedim>::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.template set_cell<dim>(dh_cell);
250
251 postprocessor->evaluate_scalar_field(
252 data.patch_values_scalar,
253 data.postprocessed_values[dataset]);
254 }
255 else
256 {
257 data.resize_system_vectors(n_components);
258
259 // at each point there is a vector valued function and
260 // its derivative...
261 if (update_flags & update_values)
262 this->dof_data[dataset]->get_function_values(
263 fe_patch_values,
265 real_part,
266 data.patch_values_system.solution_values);
267 if (update_flags & update_gradients)
268 this->dof_data[dataset]->get_function_gradients(
269 fe_patch_values,
271 real_part,
272 data.patch_values_system.solution_gradients);
273 if (update_flags & update_hessians)
274 this->dof_data[dataset]->get_function_hessians(
275 fe_patch_values,
277 real_part,
278 data.patch_values_system.solution_hessians);
279
280 if (update_flags & update_quadrature_points)
281 data.patch_values_system.evaluation_points =
282 fe_patch_values.get_quadrature_points();
283
284 const typename DoFHandler<dim,
285 spacedim>::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.template set_cell<dim>(dh_cell);
291
292 postprocessor->evaluate_vector_field(
293 data.patch_values_system,
294 data.postprocessed_values[dataset]);
295 }
296
297 for (unsigned int component = 0;
298 component < this->dof_data[dataset]->n_output_variables;
299 ++component)
300 {
301 switch (dim)
302 {
303 case 1:
304 for (unsigned int x = 0; x < n_points; ++x)
305 for (unsigned int y = 0; y < n_points; ++y)
306 my_patches[angle].data(offset + component,
307 x * n_points + y) =
308 data.postprocessed_values[dataset][x](
309 component);
310 break;
311
312 case 2:
313 for (unsigned int x = 0; x < n_points; ++x)
314 for (unsigned int y = 0; y < n_points; ++y)
315 for (unsigned int z = 0; z < n_points; ++z)
316 my_patches[angle].data(offset + component,
317 x * n_points *
318 n_points +
319 y * n_points + z) =
320 data.postprocessed_values[dataset]
321 [x * n_points + z](
322 component);
323 break;
324
325 default:
327 }
328 }
329 }
330 else if (n_components == 1)
331 {
332 this->dof_data[dataset]->get_function_values(
333 fe_patch_values,
335 real_part,
336 data.patch_values_scalar.solution_values);
337
338 switch (dim)
339 {
340 case 1:
341 for (unsigned int x = 0; x < n_points; ++x)
342 for (unsigned int y = 0; y < n_points; ++y)
343 my_patches[angle].data(offset, x * n_points + y) =
344 data.patch_values_scalar.solution_values[x];
345 break;
346
347 case 2:
348 for (unsigned int x = 0; x < n_points; ++x)
349 for (unsigned int y = 0; y < n_points; ++y)
350 for (unsigned int z = 0; z < n_points; ++z)
351 my_patches[angle].data(offset,
352 x * n_points * n_points +
353 y + z * n_points) =
354 data.patch_values_scalar
355 .solution_values[x * n_points + z];
356 break;
357
358 default:
360 }
361 }
362 else
363 // system of components
364 {
365 data.resize_system_vectors(n_components);
366 this->dof_data[dataset]->get_function_values(
367 fe_patch_values,
369 real_part,
370 data.patch_values_system.solution_values);
371
372 for (unsigned int component = 0; component < n_components;
373 ++component)
374 {
375 switch (dim)
376 {
377 case 1:
378 for (unsigned int x = 0; x < n_points; ++x)
379 for (unsigned int y = 0; y < n_points; ++y)
380 my_patches[angle].data(offset + component,
381 x * n_points + y) =
382 data.patch_values_system.solution_values[x](
383 component);
384 break;
385
386 case 2:
387 for (unsigned int x = 0; x < n_points; ++x)
388 for (unsigned int y = 0; y < n_points; ++y)
389 for (unsigned int z = 0; z < n_points; ++z)
390 my_patches[angle].data(offset + component,
391 x * n_points *
392 n_points +
393 y * n_points + z) =
394 data.patch_values_system
395 .solution_values[x * n_points + z](
396 component);
397 break;
398
399 default:
401 }
402 }
403 }
404 offset += this->dof_data[dataset]->n_output_variables;
405 }
406
407 // then do the cell data
408 for (unsigned int dataset = 0; dataset < this->cell_data.size();
409 ++dataset)
410 {
411 // we need to get at the number of the cell to which this face
412 // belongs in order to access the cell data. this is not readily
413 // available, so choose the following rather inefficient way:
414 Assert((*cell)->is_active(),
415 ExcMessage("Cell must be active for cell data"));
416 const unsigned int cell_number = std::distance(
417 this->triangulation->begin_active(),
419 *cell));
420 const double value =
421 this->cell_data[dataset]->get_cell_data_value(
422 cell_number,
424 real_part);
425 switch (dim)
426 {
427 case 1:
428 for (unsigned int x = 0; x < n_points; ++x)
429 for (unsigned int y = 0; y < n_points; ++y)
430 my_patches[angle].data(dataset + offset,
431 x * n_points + y) = value;
432 break;
433
434 case 2:
435 for (unsigned int x = 0; x < n_points; ++x)
436 for (unsigned int y = 0; y < n_points; ++y)
437 for (unsigned int z = 0; z < n_points; ++z)
438 my_patches[angle].data(dataset + offset,
439 x * n_points * n_points +
440 y * n_points + z) = value;
441 break;
442
443 default:
445 }
446 }
447 }
448 }
449}
450
451
452
453template <int dim, int spacedim>
454void
456 const unsigned int n_patches_per_circle,
457 const unsigned int nnnn_subdivisions)
458{
459 Assert(this->triangulation != nullptr,
461
462 const unsigned int n_subdivisions =
463 (nnnn_subdivisions != 0) ? nnnn_subdivisions : this->default_subdivisions;
464 Assert(n_subdivisions >= 1,
466 n_subdivisions));
467
468 this->validate_dataset_names();
469
470 unsigned int n_datasets = this->cell_data.size();
471 for (unsigned int i = 0; i < this->dof_data.size(); ++i)
472 n_datasets += this->dof_data[i]->n_output_variables;
473
475 for (unsigned int i = 0; i < this->dof_data.size(); ++i)
476 if (this->dof_data[i]->postprocessor)
477 update_flags |=
478 this->dof_data[i]->postprocessor->get_needed_update_flags();
479 // perhaps update_normal_vectors is present,
480 // which would only be useful on faces, but
481 // we may not use it here.
482 Assert(!(update_flags & update_normal_vectors),
483 ExcMessage("The update of normal vectors may not be requested for "
484 "evaluation of data on cells via DataPostprocessor."));
485
486 // first count the cells we want to
487 // create patches of and make sure
488 // there is enough memory for that
489 std::vector<cell_iterator> all_cells;
490 for (cell_iterator cell = first_cell(); cell != this->triangulation->end();
491 cell = next_cell(cell))
492 all_cells.push_back(cell);
493
494 // then also take into account that
495 // we want more than one patch to
496 // come out of every cell, as they
497 // are repeated around the axis of
498 // rotation
499 this->patches.clear();
500 this->patches.reserve(all_cells.size() * n_patches_per_circle);
501
502
503 std::vector<unsigned int> n_postprocessor_outputs(this->dof_data.size());
504 for (unsigned int dataset = 0; dataset < this->dof_data.size(); ++dataset)
505 if (this->dof_data[dataset]->postprocessor)
506 n_postprocessor_outputs[dataset] =
507 this->dof_data[dataset]->n_output_variables;
508 else
509 n_postprocessor_outputs[dataset] = 0;
510
511 Assert(!this->triangulation->is_mixed_mesh(), ExcNotImplemented());
512 const auto reference_cell = this->triangulation->get_reference_cells()[0];
514 thread_data(n_datasets,
515 n_subdivisions,
516 n_patches_per_circle,
517 n_postprocessor_outputs,
518 reference_cell.template get_default_linear_mapping<spacedim>(),
519 this->get_fes(),
520 update_flags);
521 std::vector<DataOutBase::Patch<patch_dim, patch_spacedim>> new_patches(
522 n_patches_per_circle);
523 for (unsigned int i = 0; i < new_patches.size(); ++i)
524 {
525 new_patches[i].n_subdivisions = n_subdivisions;
526 new_patches[i].reference_cell = ReferenceCells::get_hypercube<dim + 1>();
527
528 new_patches[i].data.reinit(
529 n_datasets, Utilities::fixed_power<patch_dim>(n_subdivisions + 1));
530 }
531
532 // now build the patches in parallel
534 all_cells.data(),
535 all_cells.data() + all_cells.size(),
536 [this](
537 const cell_iterator *cell,
539 &data,
540 std::vector<DataOutBase::Patch<patch_dim, patch_spacedim>> &my_patches) {
541 this->build_one_patch(cell, data, my_patches);
542 },
543 [this](const std::vector<DataOutBase::Patch<patch_dim, patch_spacedim>>
544 &new_patches) {
545 internal::DataOutRotationImplementation::append_patch_to_list<dim,
546 spacedim>(
547 new_patches, this->patches);
548 },
549 thread_data,
550 new_patches);
551}
552
553
554
555template <int dim, int spacedim>
558{
559 return this->triangulation->begin_active();
560}
561
562
563template <int dim, int spacedim>
566{
567 // convert the iterator to an
568 // active_iterator and advance
569 // this to the next active cell
571 cell;
572 ++active_cell;
573 return active_cell;
574}
575
576
577
578// explicit instantiations
579#include "numerics/data_out_rotation.inst"
580
581
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:318
Definition point.h:111
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
#define DEAL_II_NOT_IMPLEMENTED()
static ::ExceptionBase & ExcInvalidNumberOfSubdivisions(int arg1)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcNoTriangulationSelected()
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.
std::vector< index_type > data
Definition mpi.cc:734
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)
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)
constexpr double PI
Definition numbers.h:240
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices()
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)