Reference documentation for deal.II version GIT d8dacc551e 2022-08-19 06:50:03+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\}}\)
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>
23 #include <deal.II/fe/fe_values.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 
39 namespace 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 
88 template <int dim, int spacedim>
90  : surface_only(so)
91 {}
92 
93 
94 
95 template <int dim, int spacedim>
96 void
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(
211  data.patch_values_scalar,
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(
255  data.patch_values_system,
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(),
306  ExcMessage(
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 
325 template <int dim, int spacedim>
326 void
327 DataOutFaces<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 
342 template <int dim, int spacedim>
343 void
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  },
416  [this](const DataOutBase::Patch<patch_dim, patch_spacedim> &patch) {
417  internal::DataOutFacesImplementation::append_patch_to_list<dim, spacedim>(
418  patch, this->patches);
419  },
420  thread_data,
421  sample_patch);
422 }
423 
424 
425 
426 template <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 
445 template <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 void evaluate_scalar_field(const DataPostprocessorInputs::Scalar< dim > &input_data, std::vector< Vector< double >> &computed_quantities) const
virtual void evaluate_vector_field(const DataPostprocessorInputs::Vector< dim > &input_data, std::vector< Vector< double >> &computed_quantities) const
virtual UpdateFlags get_needed_update_flags() const =0
const unsigned int n_quadrature_points
Definition: fe_values.h:2432
const std::vector< Tensor< 1, spacedim > > & get_normal_vectors() const
Definition: fe_values.cc:3973
const FiniteElement< dim, spacedim > & get_fe() const
const std::vector< Point< spacedim > > & get_quadrature_points() const
unsigned int n_components() const
Definition: point.h:111
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcInvalidNumberOfSubdivisions(int arg1)
static ::ExceptionBase & ExcNoTriangulationSelected()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
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.
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:1274
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
Point< spacedim > vertices[GeometryInfo< dim >::vertices_per_cell]
static const unsigned int space_dim
unsigned int n_subdivisions
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)
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)
const FEValuesBase< dim, spacedim > & get_present_fe_values(const unsigned int dataset) const
void resize_system_vectors(const unsigned int n_components)