Reference documentation for deal.II version GIT relicensing-439-g5fda5c893d 2024-04-20 06:50: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\}}\)
Loading...
Searching...
No Matches
data_out_faces.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2000 - 2023 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
36
37
38namespace internal
39{
40 namespace DataOutFacesImplementation
41 {
42 template <int dim, int spacedim>
44 const unsigned int n_datasets,
45 const unsigned int n_subdivisions,
46 const std::vector<unsigned int> &n_postprocessor_outputs,
47 const Mapping<dim, spacedim> &mapping,
48 const std::vector<
49 std::shared_ptr<::hp::FECollection<dim, spacedim>>>
50 &finite_elements,
51 const UpdateFlags update_flags)
52 : internal::DataOutImplementation::ParallelDataBase<dim, spacedim>(
53 n_datasets,
54 n_subdivisions,
55 n_postprocessor_outputs,
56 mapping,
57 finite_elements,
58 update_flags,
59 true)
60 {}
61
62
63
68 template <int dim, int spacedim>
69 void
73 &patch,
74 std::vector<
77 &patches)
78 {
79 patches.push_back(patch);
80 patches.back().patch_index = patches.size() - 1;
81 }
82 } // namespace DataOutFacesImplementation
83} // namespace internal
84
85
86
87template <int dim, int spacedim>
89 : surface_only(so)
90{}
91
92
93
94template <int dim, int spacedim>
95void
97 const FaceDescriptor *cell_and_face,
100{
101 const cell_iterator cell = cell_and_face->first;
102 const unsigned int face_number = cell_and_face->second;
103
104 Assert(cell->is_locally_owned(), ExcNotImplemented());
105
106 // First set the kind of object we are dealing with here in the 'patch'
107 // object.
108 patch.reference_cell = cell->face(face_number)->reference_cell();
109
110 // We use the mapping to transform the vertices. However, the mapping works
111 // on cells, not faces, so transform the face vertex to a cell vertex, that
112 // to a unit cell vertex and then, finally, that to the mapped vertex. In
113 // most cases this complicated procedure will be the identity.
114 for (const unsigned int vertex : cell->face(face_number)->vertex_indices())
115 {
116 const Point<dim> vertex_reference_coordinates =
117 cell->reference_cell().template vertex<dim>(
118 cell->reference_cell().face_to_cell_vertices(
119 face_number, vertex, cell->combined_face_orientation(face_number)));
120
121 const Point<dim> vertex_real_coordinates =
122 data.mapping_collection[0].transform_unit_to_real_cell(
123 cell, vertex_reference_coordinates);
124
125 patch.vertices[vertex] = vertex_real_coordinates;
126 }
127
128
129 if (data.n_datasets > 0)
130 {
131 data.reinit_all_fe_values(this->dof_data, cell, face_number);
132 const FEValuesBase<dim> &fe_patch_values = data.get_present_fe_values(0);
133
134 const unsigned int n_q_points = fe_patch_values.n_quadrature_points;
135
136 // store the intermediate points
137 Assert(patch.space_dim == dim, ExcInternalError());
138 const std::vector<Point<dim>> &q_points =
139 fe_patch_values.get_quadrature_points();
140 // size the patch.data member in order to have enough memory for the
141 // quadrature points as well
142 patch.data.reinit(data.n_datasets + dim, q_points.size());
143 // set the flag indicating that for this cell the points are explicitly
144 // given
145 patch.points_are_available = true;
146 // copy points to patch.data
147 for (unsigned int i = 0; i < dim; ++i)
148 for (unsigned int q = 0; q < n_q_points; ++q)
149 patch.data(patch.data.size(0) - dim + i, q) = q_points[q][i];
150
151 // counter for data records
152 unsigned int offset = 0;
153
154 // first fill dof_data
155 for (unsigned int dataset = 0; dataset < this->dof_data.size(); ++dataset)
156 {
157 const FEValuesBase<dim> &this_fe_patch_values =
158 data.get_present_fe_values(dataset);
159 const unsigned int n_components =
160 this_fe_patch_values.get_fe().n_components();
161 const DataPostprocessor<dim> *postprocessor =
162 this->dof_data[dataset]->postprocessor;
163 if (postprocessor != nullptr)
164 {
165 // we have to postprocess the data, so determine, which fields
166 // have to be updated
167 const UpdateFlags update_flags =
168 postprocessor->get_needed_update_flags();
169
170 if (n_components == 1)
171 {
172 // at each point there is only one component of value,
173 // gradient etc.
174 if (update_flags & update_values)
175 this->dof_data[dataset]->get_function_values(
176 this_fe_patch_values,
178 real_part,
179 data.patch_values_scalar.solution_values);
180 if (update_flags & update_gradients)
181 this->dof_data[dataset]->get_function_gradients(
182 this_fe_patch_values,
184 real_part,
185 data.patch_values_scalar.solution_gradients);
186 if (update_flags & update_hessians)
187 this->dof_data[dataset]->get_function_hessians(
188 this_fe_patch_values,
190 real_part,
191 data.patch_values_scalar.solution_hessians);
192
193 if (update_flags & update_quadrature_points)
194 data.patch_values_scalar.evaluation_points =
195 this_fe_patch_values.get_quadrature_points();
196
197 if (update_flags & update_normal_vectors)
198 data.patch_values_scalar.normals =
199 this_fe_patch_values.get_normal_vectors();
200
202 dh_cell(&cell->get_triangulation(),
203 cell->level(),
204 cell->index(),
205 this->dof_data[dataset]->dof_handler);
206 data.patch_values_scalar.template set_cell_and_face<dim>(
207 dh_cell, face_number);
208
209 postprocessor->evaluate_scalar_field(
211 data.postprocessed_values[dataset]);
212 }
213 else
214 {
215 // at each point there is a vector valued function and its
216 // derivative...
217 data.resize_system_vectors(n_components);
218 if (update_flags & update_values)
219 this->dof_data[dataset]->get_function_values(
220 this_fe_patch_values,
222 real_part,
223 data.patch_values_system.solution_values);
224 if (update_flags & update_gradients)
225 this->dof_data[dataset]->get_function_gradients(
226 this_fe_patch_values,
228 real_part,
229 data.patch_values_system.solution_gradients);
230 if (update_flags & update_hessians)
231 this->dof_data[dataset]->get_function_hessians(
232 this_fe_patch_values,
234 real_part,
235 data.patch_values_system.solution_hessians);
236
237 if (update_flags & update_quadrature_points)
238 data.patch_values_system.evaluation_points =
239 this_fe_patch_values.get_quadrature_points();
240
241 if (update_flags & update_normal_vectors)
242 data.patch_values_system.normals =
243 this_fe_patch_values.get_normal_vectors();
244
246 dh_cell(&cell->get_triangulation(),
247 cell->level(),
248 cell->index(),
249 this->dof_data[dataset]->dof_handler);
250 data.patch_values_system.template set_cell_and_face<dim>(
251 dh_cell, face_number);
252
253 postprocessor->evaluate_vector_field(
255 data.postprocessed_values[dataset]);
256 }
257
258 for (unsigned int q = 0; q < n_q_points; ++q)
259 for (unsigned int component = 0;
260 component < this->dof_data[dataset]->n_output_variables;
261 ++component)
262 patch.data(offset + component, q) =
263 data.postprocessed_values[dataset][q](component);
264 }
265 else
266 // now we use the given data vector without modifications. again,
267 // we treat single component functions separately for efficiency
268 // reasons.
269 if (n_components == 1)
270 {
271 this->dof_data[dataset]->get_function_values(
272 this_fe_patch_values,
274 real_part,
275 data.patch_values_scalar.solution_values);
276 for (unsigned int q = 0; q < n_q_points; ++q)
277 patch.data(offset, q) =
278 data.patch_values_scalar.solution_values[q];
279 }
280 else
281 {
282 data.resize_system_vectors(n_components);
283 this->dof_data[dataset]->get_function_values(
284 this_fe_patch_values,
286 real_part,
287 data.patch_values_system.solution_values);
288 for (unsigned int component = 0; component < n_components;
289 ++component)
290 for (unsigned int q = 0; q < n_q_points; ++q)
291 patch.data(offset + component, q) =
292 data.patch_values_system.solution_values[q](component);
293 }
294 // increment the counter for the actual data record
295 offset += this->dof_data[dataset]->n_output_variables;
296 }
297
298 // then do the cell data
299 for (unsigned int dataset = 0; dataset < this->cell_data.size();
300 ++dataset)
301 {
302 // we need to get at the number of the cell to which this face
303 // belongs in order to access the cell data. this is not readily
304 // available, so choose the following rather inefficient way:
305 Assert(
306 cell->is_active(),
308 "The current function is trying to generate cell-data output "
309 "for a face that does not belong to an active cell. This is "
310 "not supported."));
311 const unsigned int cell_number = std::distance(
312 this->triangulation->begin_active(),
314
315 const double value = this->cell_data[dataset]->get_cell_data_value(
316 cell_number,
318 for (unsigned int q = 0; q < n_q_points; ++q)
319 patch.data(dataset + offset, q) = value;
320 }
321 }
322}
323
324
325
326template <int dim, int spacedim>
327void
328DataOutFaces<dim, spacedim>::build_patches(const unsigned int n_subdivisions)
329{
330 if (this->triangulation->get_reference_cells().size() == 1)
331 build_patches(this->triangulation->get_reference_cells()[0]
332 .template get_default_linear_mapping<dim, spacedim>(),
333 n_subdivisions);
334 else
335 Assert(false,
336 ExcMessage("The DataOutFaces class can currently not be "
337 "used on meshes that do not have the same cell type "
338 "throughout."));
339}
340
341
342
343template <int dim, int spacedim>
344void
346 const Mapping<dim, spacedim> &mapping,
347 const unsigned int n_subdivisions_)
348{
349 const unsigned int n_subdivisions =
350 (n_subdivisions_ != 0) ? n_subdivisions_ : this->default_subdivisions;
351
352 Assert(n_subdivisions >= 1,
354 n_subdivisions));
355
356 Assert(this->triangulation != nullptr,
358
359 this->validate_dataset_names();
360
361 unsigned int n_datasets = this->cell_data.size();
362 for (unsigned int i = 0; i < this->dof_data.size(); ++i)
363 n_datasets += this->dof_data[i]->n_output_variables;
364
365 // first collect the cells we want to create patches of; we will
366 // then iterate over them. the end-condition of the loop needs to
367 // test that next_face() returns an end iterator, as well as for the
368 // case that first_face() returns an invalid FaceDescriptor object
369 std::vector<FaceDescriptor> all_faces;
370 for (FaceDescriptor face = first_face();
371 ((face.first != this->triangulation->end()) &&
372 (face != FaceDescriptor()));
373 face = next_face(face))
374 all_faces.push_back(face);
375
376 // clear the patches array and allocate the right number of elements
377 this->patches.clear();
378 this->patches.reserve(all_faces.size());
379 Assert(this->patches.empty(), ExcInternalError());
380
381
382 std::vector<unsigned int> n_postprocessor_outputs(this->dof_data.size());
383 for (unsigned int dataset = 0; dataset < this->dof_data.size(); ++dataset)
384 if (this->dof_data[dataset]->postprocessor)
385 n_postprocessor_outputs[dataset] =
386 this->dof_data[dataset]->n_output_variables;
387 else
388 n_postprocessor_outputs[dataset] = 0;
389
390 UpdateFlags update_flags = update_values;
391 for (unsigned int i = 0; i < this->dof_data.size(); ++i)
392 if (this->dof_data[i]->postprocessor)
393 update_flags |=
394 this->dof_data[i]->postprocessor->get_needed_update_flags();
395 update_flags |= update_quadrature_points;
396
398 n_datasets,
399 n_subdivisions,
400 n_postprocessor_outputs,
401 mapping,
402 this->get_fes(),
403 update_flags);
405 sample_patch.n_subdivisions = n_subdivisions;
406
407 // now build the patches in parallel
409 all_faces.data(),
410 all_faces.data() + all_faces.size(),
411 [this](
412 const FaceDescriptor *cell_and_face,
415 this->build_one_patch(cell_and_face, data, patch);
416 },
418 internal::DataOutFacesImplementation::append_patch_to_list<dim, spacedim>(
419 patch, this->patches);
420 },
421 thread_data,
422 sample_patch);
423}
424
425
426
427template <int dim, int spacedim>
430{
431 // simply find first active cell with a face on the boundary
432 for (const auto &cell : this->triangulation->active_cell_iterators())
433 if (cell->is_locally_owned())
434 for (const unsigned int f : cell->face_indices())
435 if ((surface_only == false) || cell->face(f)->at_boundary())
436 return FaceDescriptor(cell, f);
437
438 // just return an invalid descriptor if we haven't found a locally
439 // owned face. this can happen in parallel where all boundary
440 // faces are owned by other processors
441 return FaceDescriptor();
442}
443
444
445
446template <int dim, int spacedim>
449{
450 FaceDescriptor face = old_face;
451
452 // first check whether the present cell has more faces on the boundary. since
453 // we started with this face, its cell must clearly be locally owned
454 Assert(face.first->is_locally_owned(), ExcInternalError());
455 for (unsigned int f = face.second + 1; f < face.first->n_faces(); ++f)
456 if (!surface_only || face.first->face(f)->at_boundary())
457 // yup, that is so, so return it
458 {
459 face.second = f;
460 return face;
461 }
462
463 // otherwise find the next active cell that has a face on the boundary
464
465 // convert the iterator to an active_iterator and advance this to the next
466 // active cell
468 face.first;
469
470 // increase face pointer by one
471 ++active_cell;
472
473 // while there are active cells
474 while (active_cell != this->triangulation->end())
475 {
476 // check all the faces of this active cell. but skip it altogether
477 // if it isn't locally owned
478 if (active_cell->is_locally_owned())
479 for (const unsigned int f : face.first->face_indices())
480 if (!surface_only || active_cell->face(f)->at_boundary())
481 {
482 face.first = active_cell;
483 face.second = f;
484 return face;
485 }
486
487 // the present cell had no faces on the boundary (or was not locally
488 // owned), so check next cell
489 ++active_cell;
490 }
491
492 // we fell off the edge, so return with invalid pointer
493 face.first = this->triangulation->end();
494 face.second = 0;
495 return face;
496}
497
498
499
500// explicit instantiations
501#include "data_out_faces.inst"
502
void build_one_patch(const FaceDescriptor *cell_and_face, internal::DataOutFacesImplementation::ParallelData< dim, spacedim > &data, DataOutBase::Patch< patch_dim, patch_spacedim > &patch)
typename DataOut_DoFData< dim, patch_dim, spacedim, patch_spacedim >::cell_iterator cell_iterator
virtual FaceDescriptor next_face(const FaceDescriptor &face)
typename std::pair< cell_iterator, unsigned int > FaceDescriptor
virtual void build_patches(const unsigned int n_subdivisions=0)
DataOutFaces(const bool surface_only=true)
virtual FaceDescriptor first_face()
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 unsigned int n_quadrature_points
const std::vector< Tensor< 1, spacedim > > & get_normal_vectors() 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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static ::ExceptionBase & ExcInvalidNumberOfSubdivisions(int arg1)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNoTriangulationSelected()
static ::ExceptionBase & ExcMessage(std::string arg1)
typename ActiveSelector::active_cell_iterator active_cell_iterator
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 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 DataOutBase::Patch< DataOutFaces< dim, spacedim >::patch_dim, DataOutFaces< dim, spacedim >::patch_spacedim > &patch, std::vector< DataOutBase::Patch< DataOutFaces< dim, spacedim >::patch_dim, DataOutFaces< dim, spacedim >::patch_spacedim > > &patches)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
ReferenceCell reference_cell
Table< 2, float > data
static const unsigned int space_dim
unsigned int n_subdivisions
std::array< Point< spacedim >, GeometryInfo< dim >::vertices_per_cell > vertices
ParallelData(const unsigned int n_datasets, const unsigned int n_subdivisions, 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)
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
const ::hp::MappingCollection< dim, spacedim > mapping_collection
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)