deal.II version GIT relicensing-2167-g9622207b8f 2024-11-21 12:40:00+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: LGPL-2.1-or-later
4// Copyright (C) 2000 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
17
20
21#include <deal.II/fe/fe.h>
23#include <deal.II/fe/mapping.h>
24
25#include <deal.II/grid/tria.h>
27
29
31#include <deal.II/lac/vector.h>
32
34
35#include <cmath>
36
38
39
40// TODO: Update documentation
41// TODO: Unify code for dimensions
42
43
44// TODO: build_some_patches isn't going to work if first_cell/next_cell
45// don't iterate over all cells and if cell data is requested. in that
46// case, we need to calculate cell_number as in the DataOut class
47
48// Not implemented for 3d
49
50
51namespace internal
52{
53 namespace DataOutRotationImplementation
54 {
55 template <int dim, int spacedim>
57 const unsigned int n_datasets,
58 const unsigned int n_subdivisions,
59 const unsigned int n_patches_per_circle,
60 const std::vector<unsigned int> &n_postprocessor_outputs,
61 const Mapping<dim, spacedim> &mapping,
62 const std::vector<
63 std::shared_ptr<::hp::FECollection<dim, spacedim>>>
64 &finite_elements,
65 const UpdateFlags update_flags)
66 : internal::DataOutImplementation::ParallelDataBase<dim, spacedim>(
67 n_datasets,
68 n_subdivisions,
69 n_postprocessor_outputs,
70 mapping,
71 finite_elements,
72 update_flags,
73 false)
74 , n_patches_per_circle(n_patches_per_circle)
75 {}
76
77
78
83 template <int dim, int spacedim>
84 void
86 const std::vector<
89 &new_patches,
90 std::vector<
93 &patches)
94 {
95 for (unsigned int i = 0; i < new_patches.size(); ++i)
96 {
97 patches.push_back(new_patches[i]);
98 patches.back().patch_index = patches.size() - 1;
99 }
100 }
101 } // namespace DataOutRotationImplementation
102} // namespace internal
103
104
105
106template <int dim, int spacedim>
107void
109 const cell_iterator *cell,
111 std::vector<DataOutBase::Patch<patch_dim, patch_spacedim>> &my_patches)
112{
113 if (dim == 3)
114 {
115 // would this function make any sense after all? who would want to
116 // output/compute in four space dimensions?
118 return;
119 }
120
121 Assert((*cell)->is_locally_owned(), ExcNotImplemented());
122
123 const unsigned int n_patches_per_circle = data.n_patches_per_circle;
124
125 // another abbreviation denoting the number of q_points in each direction
126 const unsigned int n_points = data.n_subdivisions + 1;
127
128 // set up an array that holds the directions in the plane of rotation in
129 // which we will put points in the whole domain (not the rotationally
130 // reduced one in which the computation took place. for simplicity add the
131 // initial direction at the end again
132 std::vector<Point<dim + 1>> angle_directions(n_patches_per_circle + 1);
133 for (unsigned int i = 0; i <= n_patches_per_circle; ++i)
134 {
135 angle_directions[i][dim - 1] =
136 std::cos(2 * numbers::PI * i / n_patches_per_circle);
137 angle_directions[i][dim] =
138 std::sin(2 * numbers::PI * i / n_patches_per_circle);
139 }
140
141 for (unsigned int angle = 0; angle < n_patches_per_circle; ++angle)
142 {
143 // first compute the vertices of the patch. note that they will have to
144 // be computed from the vertices of the cell, which has one dim
145 // less, however.
146 switch (dim)
147 {
148 case 1:
149 {
150 const double r1 = (*cell)->vertex(0)[0],
151 r2 = (*cell)->vertex(1)[0];
152 Assert(r1 >= 0, ExcRadialVariableHasNegativeValues(r1));
153 Assert(r2 >= 0, ExcRadialVariableHasNegativeValues(r2));
154
155 my_patches[angle].vertices[0] = r1 * angle_directions[angle];
156 my_patches[angle].vertices[1] = r2 * angle_directions[angle];
157 my_patches[angle].vertices[2] = r1 * angle_directions[angle + 1];
158 my_patches[angle].vertices[3] = r2 * angle_directions[angle + 1];
159
160 break;
161 }
162
163 case 2:
164 {
165 for (const unsigned int vertex :
167 {
168 const Point<dim> v = (*cell)->vertex(vertex);
169
170 // make sure that the radial variable is nonnegative
171 Assert(v[0] >= 0, ExcRadialVariableHasNegativeValues(v[0]));
172
173 // now set the vertices of the patch
174 my_patches[angle].vertices[vertex] =
175 v[0] * angle_directions[angle];
176 my_patches[angle].vertices[vertex][0] = v[1];
177
178 my_patches[angle]
179 .vertices[vertex + GeometryInfo<dim>::vertices_per_cell] =
180 v[0] * angle_directions[angle + 1];
181 my_patches[angle]
182 .vertices[vertex + GeometryInfo<dim>::vertices_per_cell]
183 [0] = v[1];
184 }
185
186 break;
187 }
188
189 default:
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<dim> &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,
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,
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,
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 DoFHandler<dim,
246 spacedim>::active_cell_iterator
247 dh_cell(&(*cell)->get_triangulation(),
248 (*cell)->level(),
249 (*cell)->index(),
250 this->dof_data[dataset]->dof_handler);
251 data.patch_values_scalar.template set_cell<dim>(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,
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,
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,
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 DoFHandler<dim,
287 spacedim>::active_cell_iterator
288 dh_cell(&(*cell)->get_triangulation(),
289 (*cell)->level(),
290 (*cell)->index(),
291 this->dof_data[dataset]->dof_handler);
292 data.patch_values_system.template set_cell<dim>(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 (dim)
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:
329 }
330 }
331 }
332 else if (n_components == 1)
333 {
334 this->dof_data[dataset]->get_function_values(
335 fe_patch_values,
337 real_part,
338 data.patch_values_scalar.solution_values);
339
340 switch (dim)
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:
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,
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 (dim)
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:
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)->is_active(),
417 ExcMessage("Cell must be active for cell data"));
418 const unsigned int cell_number = std::distance(
419 this->triangulation->begin_active(),
421 *cell));
422 const double value =
423 this->cell_data[dataset]->get_cell_data_value(
424 cell_number,
426 real_part);
427 switch (dim)
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:
447 }
448 }
449 }
450 }
451}
452
453
454
455template <int dim, int spacedim>
456void
458 const unsigned int n_patches_per_circle,
459 const unsigned int nnnn_subdivisions)
460{
461 Assert(this->triangulation != nullptr,
463
464 const unsigned int n_subdivisions =
465 (nnnn_subdivisions != 0) ? nnnn_subdivisions : this->default_subdivisions;
466 Assert(n_subdivisions >= 1,
468 n_subdivisions));
469
470 this->validate_dataset_names();
471
472 unsigned int n_datasets = this->cell_data.size();
473 for (unsigned int i = 0; i < this->dof_data.size(); ++i)
474 n_datasets += this->dof_data[i]->n_output_variables;
475
477 for (unsigned int i = 0; i < this->dof_data.size(); ++i)
478 if (this->dof_data[i]->postprocessor)
479 update_flags |=
480 this->dof_data[i]->postprocessor->get_needed_update_flags();
481 // perhaps update_normal_vectors is present,
482 // which would only be useful on faces, but
483 // we may not use it here.
484 Assert(!(update_flags & update_normal_vectors),
485 ExcMessage("The update of normal vectors may not be requested for "
486 "evaluation of data on cells via DataPostprocessor."));
487
488 // first count the cells we want to
489 // create patches of and make sure
490 // there is enough memory for that
491 std::vector<cell_iterator> all_cells;
492 for (cell_iterator cell = first_cell(); cell != this->triangulation->end();
493 cell = next_cell(cell))
494 all_cells.push_back(cell);
495
496 // then also take into account that
497 // we want more than one patch to
498 // come out of every cell, as they
499 // are repeated around the axis of
500 // rotation
501 this->patches.clear();
502 this->patches.reserve(all_cells.size() * n_patches_per_circle);
503
504
505 std::vector<unsigned int> n_postprocessor_outputs(this->dof_data.size());
506 for (unsigned int dataset = 0; dataset < this->dof_data.size(); ++dataset)
507 if (this->dof_data[dataset]->postprocessor)
508 n_postprocessor_outputs[dataset] =
509 this->dof_data[dataset]->n_output_variables;
510 else
511 n_postprocessor_outputs[dataset] = 0;
512
513 Assert(!this->triangulation->is_mixed_mesh(), ExcNotImplemented());
514 const auto reference_cell = this->triangulation->get_reference_cells()[0];
516 thread_data(
517 n_datasets,
518 n_subdivisions,
519 n_patches_per_circle,
520 n_postprocessor_outputs,
521 reference_cell.template get_default_linear_mapping<dim, spacedim>(),
522 this->get_fes(),
523 update_flags);
524 std::vector<DataOutBase::Patch<patch_dim, patch_spacedim>> new_patches(
525 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].reference_cell = ReferenceCells::get_hypercube<dim + 1>();
530
531 new_patches[i].data.reinit(
532 n_datasets, Utilities::fixed_power<patch_dim>(n_subdivisions + 1));
533 }
534
535 // now build the patches in parallel
537 all_cells.data(),
538 all_cells.data() + all_cells.size(),
539 [this](
540 const cell_iterator *cell,
542 &data,
543 std::vector<DataOutBase::Patch<patch_dim, patch_spacedim>> &my_patches) {
544 this->build_one_patch(cell, data, my_patches);
545 },
546 [this](const std::vector<DataOutBase::Patch<patch_dim, patch_spacedim>>
547 &new_patches) {
548 internal::DataOutRotationImplementation::append_patch_to_list<dim,
549 spacedim>(
550 new_patches, this->patches);
551 },
552 thread_data,
553 new_patches);
554}
555
556
557
558template <int dim, int spacedim>
561{
562 return this->triangulation->begin_active();
563}
564
565
566template <int dim, int spacedim>
569{
570 // convert the iterator to an
571 // active_iterator and advance
572 // this to the next active cell
574 cell;
575 ++active_cell;
576 return active_cell;
577}
578
579
580
581// explicit instantiations
582#include "data_out_rotation.inst"
583
584
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:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
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.
#define DEAL_II_NOT_IMPLEMENTED()
std::vector< index_type > data
Definition mpi.cc:735
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)
static constexpr double PI
Definition numbers.h:254
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
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)