Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
data_out_faces.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 2019 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 
16 #include <deal.II/base/quadrature_lib.h>
17 #include <deal.II/base/work_stream.h>
18 
19 #include <deal.II/dofs/dof_accessor.h>
20 #include <deal.II/dofs/dof_handler.h>
21 
22 #include <deal.II/fe/fe.h>
23 #include <deal.II/fe/fe_values.h>
24 #include <deal.II/fe/mapping_q1.h>
25 
26 #include <deal.II/grid/tria.h>
27 #include <deal.II/grid/tria_iterator.h>
28 
29 #include <deal.II/hp/fe_values.h>
30 
31 #include <deal.II/lac/block_vector.h>
32 #include <deal.II/lac/vector.h>
33 
34 #include <deal.II/numerics/data_out_faces.h>
35 
36 DEAL_II_NAMESPACE_OPEN
37 
38 
39 namespace internal
40 {
41  namespace DataOutFacesImplementation
42  {
43  template <int dim, int spacedim>
44  ParallelData<dim, spacedim>::ParallelData(
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
71  append_patch_to_list(
73  std::vector<DataOutBase::Patch<dim - 1, spacedim>> &patches)
74  {
75  patches.push_back(patch);
76  patches.back().patch_index = patches.size() - 1;
77  }
78  } // namespace DataOutFacesImplementation
79 } // namespace internal
80 
81 
82 
83 template <int dim, typename DoFHandlerType>
85  : surface_only(so)
86 {
87  Assert(dim == DoFHandlerType::dimension, ExcNotImplemented());
88 }
89 
90 
91 
92 template <int dim, typename DoFHandlerType>
93 void
95  const FaceDescriptor *cell_and_face,
97  & data,
99 {
100  Assert(cell_and_face->first->is_locally_owned(), ExcNotImplemented());
101 
102  // we use the mapping to transform the vertices. However, the mapping works
103  // on cells, not faces, so transform the face vertex to a cell vertex, that
104  // to a unit cell vertex and then, finally, that to the mapped vertex. In
105  // most cases this complicated procedure will be the identity.
106  for (unsigned int vertex = 0;
107  vertex < GeometryInfo<dimension - 1>::vertices_per_cell;
108  ++vertex)
109  patch.vertices[vertex] =
110  data.mapping_collection[0].transform_unit_to_real_cell(
111  cell_and_face->first,
114  cell_and_face->second,
115  vertex,
116  cell_and_face->first->face_orientation(cell_and_face->second),
117  cell_and_face->first->face_flip(cell_and_face->second),
118  cell_and_face->first->face_rotation(cell_and_face->second))));
119 
120  if (data.n_datasets > 0)
121  {
122  data.reinit_all_fe_values(this->dof_data,
123  cell_and_face->first,
124  cell_and_face->second);
125  const FEValuesBase<dimension> &fe_patch_values =
126  data.get_present_fe_values(0);
127 
128  const unsigned int n_q_points = fe_patch_values.n_quadrature_points;
129 
130  // store the intermediate points
131  Assert(patch.space_dim == dimension, ExcInternalError());
132  const std::vector<Point<dimension>> &q_points =
133  fe_patch_values.get_quadrature_points();
134  // resize the patch.data member in order to have enough memory for the
135  // quadrature points as well
136  patch.data.reinit(data.n_datasets + dimension, patch.data.size(1));
137  // set the flag indicating that for this cell the points are explicitly
138  // given
139  patch.points_are_available = true;
140  // copy points to patch.data
141  for (unsigned int i = 0; i < dimension; ++i)
142  for (unsigned int q = 0; q < n_q_points; ++q)
143  patch.data(patch.data.size(0) - dimension + i, q) = q_points[q][i];
144 
145  // counter for data records
146  unsigned int offset = 0;
147 
148  // first fill dof_data
149  for (unsigned int dataset = 0; dataset < this->dof_data.size(); ++dataset)
150  {
151  const FEValuesBase<dimension> &this_fe_patch_values =
152  data.get_present_fe_values(dataset);
153  const unsigned int n_components =
154  this_fe_patch_values.get_fe().n_components();
155  const DataPostprocessor<dim> *postprocessor =
156  this->dof_data[dataset]->postprocessor;
157  if (postprocessor != nullptr)
158  {
159  // we have to postprocess the data, so determine, which fields
160  // have to be updated
161  const UpdateFlags update_flags =
162  postprocessor->get_needed_update_flags();
163 
164  if (n_components == 1)
165  {
166  // at each point there is only one component of value,
167  // gradient etc.
168  if (update_flags & update_values)
169  this->dof_data[dataset]->get_function_values(
170  this_fe_patch_values,
171  internal::DataOutImplementation::ComponentExtractor::
172  real_part,
173  data.patch_values_scalar.solution_values);
174  if (update_flags & update_gradients)
175  this->dof_data[dataset]->get_function_gradients(
176  this_fe_patch_values,
177  internal::DataOutImplementation::ComponentExtractor::
178  real_part,
179  data.patch_values_scalar.solution_gradients);
180  if (update_flags & update_hessians)
181  this->dof_data[dataset]->get_function_hessians(
182  this_fe_patch_values,
183  internal::DataOutImplementation::ComponentExtractor::
184  real_part,
185  data.patch_values_scalar.solution_hessians);
186 
187  if (update_flags & update_quadrature_points)
188  data.patch_values_scalar.evaluation_points =
189  this_fe_patch_values.get_quadrature_points();
190 
191  if (update_flags & update_normal_vectors)
192  data.patch_values_scalar.normals =
193  this_fe_patch_values.get_all_normal_vectors();
194 
195  const typename DoFHandlerType::active_cell_iterator dh_cell(
196  &cell_and_face->first->get_triangulation(),
197  cell_and_face->first->level(),
198  cell_and_face->first->index(),
199  this->dof_data[dataset]->dof_handler);
200  data.patch_values_scalar.template set_cell<DoFHandlerType>(
201  dh_cell);
202 
203  postprocessor->evaluate_scalar_field(
204  data.patch_values_scalar,
205  data.postprocessed_values[dataset]);
206  }
207  else
208  {
209  // at each point there is a vector valued function and its
210  // derivative...
211  data.resize_system_vectors(n_components);
212  if (update_flags & update_values)
213  this->dof_data[dataset]->get_function_values(
214  this_fe_patch_values,
215  internal::DataOutImplementation::ComponentExtractor::
216  real_part,
217  data.patch_values_system.solution_values);
218  if (update_flags & update_gradients)
219  this->dof_data[dataset]->get_function_gradients(
220  this_fe_patch_values,
221  internal::DataOutImplementation::ComponentExtractor::
222  real_part,
223  data.patch_values_system.solution_gradients);
224  if (update_flags & update_hessians)
225  this->dof_data[dataset]->get_function_hessians(
226  this_fe_patch_values,
227  internal::DataOutImplementation::ComponentExtractor::
228  real_part,
229  data.patch_values_system.solution_hessians);
230 
231  if (update_flags & update_quadrature_points)
232  data.patch_values_system.evaluation_points =
233  this_fe_patch_values.get_quadrature_points();
234 
235  if (update_flags & update_normal_vectors)
236  data.patch_values_system.normals =
237  this_fe_patch_values.get_all_normal_vectors();
238 
239  const typename DoFHandlerType::active_cell_iterator dh_cell(
240  &cell_and_face->first->get_triangulation(),
241  cell_and_face->first->level(),
242  cell_and_face->first->index(),
243  this->dof_data[dataset]->dof_handler);
244  data.patch_values_system.template set_cell<DoFHandlerType>(
245  dh_cell);
246 
247  postprocessor->evaluate_vector_field(
248  data.patch_values_system,
249  data.postprocessed_values[dataset]);
250  }
251 
252  for (unsigned int q = 0; q < n_q_points; ++q)
253  for (unsigned int component = 0;
254  component < this->dof_data[dataset]->n_output_variables;
255  ++component)
256  patch.data(offset + component, q) =
257  data.postprocessed_values[dataset][q](component);
258  }
259  else
260  // now we use the given data vector without modifications. again,
261  // we treat single component functions separately for efficiency
262  // reasons.
263  if (n_components == 1)
264  {
265  this->dof_data[dataset]->get_function_values(
266  this_fe_patch_values,
267  internal::DataOutImplementation::ComponentExtractor::real_part,
268  data.patch_values_scalar.solution_values);
269  for (unsigned int q = 0; q < n_q_points; ++q)
270  patch.data(offset, q) =
271  data.patch_values_scalar.solution_values[q];
272  }
273  else
274  {
275  data.resize_system_vectors(n_components);
276  this->dof_data[dataset]->get_function_values(
277  this_fe_patch_values,
278  internal::DataOutImplementation::ComponentExtractor::real_part,
279  data.patch_values_system.solution_values);
280  for (unsigned int component = 0; component < n_components;
281  ++component)
282  for (unsigned int q = 0; q < n_q_points; ++q)
283  patch.data(offset + component, q) =
284  data.patch_values_system.solution_values[q](component);
285  }
286  // increment the counter for the actual data record
287  offset += this->dof_data[dataset]->n_output_variables;
288  }
289 
290  // then do the cell data
291  for (unsigned int dataset = 0; dataset < this->cell_data.size();
292  ++dataset)
293  {
294  // we need to get at the number of the cell to which this face
295  // belongs in order to access the cell data. this is not readily
296  // available, so choose the following rather inefficient way:
297  Assert(
298  cell_and_face->first->active(),
299  ExcMessage(
300  "The current function is trying to generate cell-data output "
301  "for a face that does not belong to an active cell. This is "
302  "not supported."));
303  const unsigned int cell_number =
304  std::distance(this->triangulation->begin_active(),
306  active_cell_iterator(cell_and_face->first));
307 
308  const double value = this->cell_data[dataset]->get_cell_data_value(
309  cell_number,
310  internal::DataOutImplementation::ComponentExtractor::real_part);
311  for (unsigned int q = 0; q < n_q_points; ++q)
312  patch.data(dataset + offset, q) = value;
313  }
314  }
315 }
316 
317 
318 
319 template <int dim, typename DoFHandlerType>
320 void
322  const unsigned int n_subdivisions_)
323 {
324  build_patches(StaticMappingQ1<dimension>::mapping, n_subdivisions_);
325 }
326 
327 
328 
329 template <int dim, typename DoFHandlerType>
330 void
332  const Mapping<dimension> &mapping,
333  const unsigned int n_subdivisions_)
334 {
335  // Check consistency of redundant template parameter
336  Assert(dim == dimension, ExcDimensionMismatch(dim, dimension));
337 
338  const unsigned int n_subdivisions =
339  (n_subdivisions_ != 0) ? n_subdivisions_ : this->default_subdivisions;
340 
341  Assert(n_subdivisions >= 1,
343  n_subdivisions));
344 
345  Assert(this->triangulation != nullptr,
347 
348  this->validate_dataset_names();
349 
350  unsigned int n_datasets = this->cell_data.size();
351  for (unsigned int i = 0; i < this->dof_data.size(); ++i)
352  n_datasets += this->dof_data[i]->n_output_variables;
353 
354  // first collect the cells we want to create patches of; we will
355  // then iterate over them. the end-condition of the loop needs to
356  // test that next_face() returns an end iterator, as well as for the
357  // case that first_face() returns an invalid FaceDescriptor object
358  std::vector<FaceDescriptor> all_faces;
359  for (FaceDescriptor face = first_face();
360  ((face.first != this->triangulation->end()) &&
361  (face != FaceDescriptor()));
362  face = next_face(face))
363  all_faces.push_back(face);
364 
365  // clear the patches array and allocate the right number of elements
366  this->patches.clear();
367  this->patches.reserve(all_faces.size());
368  Assert(this->patches.size() == 0, ExcInternalError());
369 
370 
371  std::vector<unsigned int> n_postprocessor_outputs(this->dof_data.size());
372  for (unsigned int dataset = 0; dataset < this->dof_data.size(); ++dataset)
373  if (this->dof_data[dataset]->postprocessor)
374  n_postprocessor_outputs[dataset] =
375  this->dof_data[dataset]->n_output_variables;
376  else
377  n_postprocessor_outputs[dataset] = 0;
378 
379  UpdateFlags update_flags = update_values;
380  for (unsigned int i = 0; i < this->dof_data.size(); ++i)
381  if (this->dof_data[i]->postprocessor)
382  update_flags |=
383  this->dof_data[i]->postprocessor->get_needed_update_flags();
384  update_flags |= update_quadrature_points;
385 
387  thread_data(n_datasets,
388  n_subdivisions,
389  n_postprocessor_outputs,
390  mapping,
391  this->get_fes(),
392  update_flags);
393  DataOutBase::Patch<dimension - 1, space_dimension> sample_patch;
394  sample_patch.n_subdivisions = n_subdivisions;
395  sample_patch.data.reinit(
396  n_datasets, Utilities::fixed_power<dimension - 1>(n_subdivisions + 1));
397 
398  // now build the patches in parallel
399  WorkStream::run(all_faces.data(),
400  all_faces.data() + all_faces.size(),
402  this,
403  std::placeholders::_1,
404  std::placeholders::_2,
405  std::placeholders::_3),
407  append_patch_to_list<dim, space_dimension>,
408  std::placeholders::_1,
409  std::ref(this->patches)),
410  thread_data,
411  sample_patch);
412 }
413 
414 
415 
416 template <int dim, typename DoFHandlerType>
419 {
420  // simply find first active cell with a face on the boundary
422  cell = this->triangulation->begin_active();
423  for (; cell != this->triangulation->end(); ++cell)
424  if (cell->is_locally_owned())
425  for (unsigned int f = 0; f < GeometryInfo<dimension>::faces_per_cell; ++f)
426  if (!surface_only || cell->face(f)->at_boundary())
427  return FaceDescriptor(cell, f);
428 
429  // just return an invalid descriptor if we haven't found a locally
430  // owned face. this can happen in parallel where all boundary
431  // faces are owned by other processors
432  return FaceDescriptor();
433 }
434 
435 
436 
437 template <int dim, typename DoFHandlerType>
440 {
441  FaceDescriptor face = old_face;
442 
443  // first check whether the present cell has more faces on the boundary. since
444  // we started with this face, its cell must clearly be locally owned
445  Assert(face.first->is_locally_owned(), ExcInternalError());
446  for (unsigned int f = face.second + 1;
448  ++f)
449  if (!surface_only || face.first->face(f)->at_boundary())
450  // yup, that is so, so return it
451  {
452  face.second = f;
453  return face;
454  }
455 
456  // otherwise find the next active cell that has a face on the boundary
457 
458  // convert the iterator to an active_iterator and advance this to the next
459  // active cell
461  active_cell = face.first;
462 
463  // increase face pointer by one
464  ++active_cell;
465 
466  // while there are active cells
467  while (active_cell != this->triangulation->end())
468  {
469  // check all the faces of this active cell. but skip it altogether
470  // if it isn't locally owned
471  if (active_cell->is_locally_owned())
472  for (unsigned int f = 0; f < GeometryInfo<dimension>::faces_per_cell;
473  ++f)
474  if (!surface_only || active_cell->face(f)->at_boundary())
475  {
476  face.first = active_cell;
477  face.second = f;
478  return face;
479  }
480 
481  // the present cell had no faces on the boundary (or was not locally
482  // owned), so check next cell
483  ++active_cell;
484  }
485 
486  // we fell off the edge, so return with invalid pointer
487  face.first = this->triangulation->end();
488  face.second = 0;
489  return face;
490 }
491 
492 
493 
494 // explicit instantiations
495 #include "data_out_faces.inst"
496 
497 DEAL_II_NAMESPACE_CLOSE
Shape function values.
static ::ExceptionBase & ExcNoTriangulationSelected()
const FiniteElement< dim, spacedim > & get_fe() const
virtual void evaluate_vector_field(const DataPostprocessorInputs::Vector< dim > &input_data, std::vector< Vector< double >> &computed_quantities) const
Transformed quadrature points.
TriaActiveIterator< CellAccessor< dim, spacedim > > active_cell_iterator
Definition: tria.h:1529
static ::ExceptionBase & ExcInvalidNumberOfSubdivisions(int arg1)
Point< spacedim > vertices[GeometryInfo< dim >::vertices_per_cell]
virtual FaceDescriptor first_face()
const std::vector< Point< spacedim > > & get_quadrature_points() const
static ::ExceptionBase & ExcMessage(std::string arg1)
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
#define Assert(cond, exc)
Definition: exceptions.h:1407
UpdateFlags
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Abstract base class for mapping classes.
Definition: dof_tools.h:57
virtual FaceDescriptor next_face(const FaceDescriptor &face)
unsigned int n_subdivisions
Second derivatives of shape functions.
const std::vector< Tensor< 1, spacedim > > & get_all_normal_vectors() const
Definition: fe_values.cc:4205
DataOutFaces(const bool surface_only=true)
const unsigned int n_quadrature_points
Definition: fe_values.h:2099
size_type size(const unsigned int i) const
static const unsigned int space_dim
Shape function gradients.
virtual void evaluate_scalar_field(const DataPostprocessorInputs::Scalar< dim > &input_data, std::vector< Vector< double >> &computed_quantities) const
Normal vectors.
typename std::pair< cell_iterator, unsigned int > FaceDescriptor
static ::ExceptionBase & ExcNotImplemented()
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:1167
void build_one_patch(const FaceDescriptor *cell_and_face, internal::DataOutFacesImplementation::ParallelData< dimension, dimension > &data, DataOutBase::Patch< dimension - 1, space_dimension > &patch)
virtual void build_patches(const unsigned int n_subdivisions=0)
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
Table< 2, float > data
static ::ExceptionBase & ExcInternalError()
virtual UpdateFlags get_needed_update_flags() const =0