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