Reference documentation for deal.II version 8.4.1
data_out.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2015 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 at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/work_stream.h>
17 #include <deal.II/numerics/data_out.h>
18 #include <deal.II/grid/tria.h>
19 #include <deal.II/dofs/dof_handler.h>
20 #include <deal.II/dofs/dof_accessor.h>
21 #include <deal.II/hp/dof_handler.h>
22 #include <deal.II/grid/tria_iterator.h>
23 #include <deal.II/fe/fe.h>
24 #include <deal.II/fe/fe_dgq.h>
25 #include <deal.II/fe/fe_values.h>
26 #include <deal.II/hp/fe_values.h>
27 #include <deal.II/fe/mapping_q1.h>
28 
29 #include <sstream>
30 
31 DEAL_II_NAMESPACE_OPEN
32 
33 
34 namespace internal
35 {
36  namespace DataOut
37  {
38  template <int dim, int spacedim>
39  ParallelData<dim,spacedim>::
40  ParallelData (const unsigned int n_datasets,
41  const unsigned int n_subdivisions,
42  const std::vector<unsigned int> &n_postprocessor_outputs,
43  const Mapping<dim,spacedim> &mapping,
44  const std::vector<std_cxx11::shared_ptr<::hp::FECollection<dim,spacedim> > > &finite_elements,
45  const UpdateFlags update_flags,
46  const std::vector<std::vector<unsigned int> > &cell_to_patch_index_map)
47  :
48  ParallelDataBase<dim,spacedim> (n_datasets,
49  n_subdivisions,
50  n_postprocessor_outputs,
51  mapping,
52  finite_elements,
53  update_flags,
54  false),
55  cell_to_patch_index_map (&cell_to_patch_index_map)
56  {}
57  }
58 }
59 
60 
61 
62 template <int dim, typename DoFHandlerType>
63 void
66 (const std::pair<cell_iterator, unsigned int> *cell_and_index,
68  const unsigned int n_subdivisions,
69  const CurvedCellRegion curved_cell_region,
71 {
72  // first create the output object that we will write into
74  patch.n_subdivisions = n_subdivisions;
75 
76  // use ucd_to_deal map as patch vertices are in the old, unnatural
77  // ordering. if the mapping does not preserve locations
78  // (e.g. MappingQEulerian), we need to compute the offset of the vertex for
79  // the graphical output. Otherwise, we can just use the vertex info.
80  for (unsigned int vertex=0; vertex<GeometryInfo<DoFHandlerType::dimension>::vertices_per_cell; ++vertex)
81  if (scratch_data.mapping_collection[0].preserves_vertex_locations())
82  patch.vertices[vertex] = cell_and_index->first->vertex(vertex);
83  else
84  patch.vertices[vertex] = scratch_data.mapping_collection[0].transform_unit_to_real_cell
85  (cell_and_index->first,
87 
88  if (scratch_data.n_datasets > 0)
89  {
90  // create DoFHandlerType::active_cell_iterator and initialize FEValues
91  scratch_data.reinit_all_fe_values(this->dof_data, cell_and_index->first);
92 
94  = scratch_data.get_present_fe_values (0);
95 
96  const unsigned int n_q_points = fe_patch_values.n_quadrature_points;
97 
98  // depending on the requested output of curved cells, if necessary
99  // append the quadrature points to the last rows of the patch.data
100  // member. This is the case if we want to produce curved cells at the
101  // boundary and this cell actually is at the boundary, or else if we
102  // want to produce curved cells everywhere
103  //
104  // note: a cell is *always* at the boundary if dim<spacedim
105  if (curved_cell_region==curved_inner_cells
106  ||
107  (curved_cell_region==curved_boundary
108  &&
109  (cell_and_index->first->at_boundary()
110  ||
111  (DoFHandlerType::dimension != DoFHandlerType::space_dimension))))
112  {
113  Assert(patch.space_dim==DoFHandlerType::space_dimension, ExcInternalError());
114  const std::vector<Point<DoFHandlerType::space_dimension> > &q_points=fe_patch_values.get_quadrature_points();
115  // resize the patch.data member in order to have enough memory for
116  // the quadrature points as well
117  patch.data.reinit (scratch_data.n_datasets+DoFHandlerType::space_dimension, n_q_points);
118  // set the flag indicating that for this cell the points are
119  // explicitly given
120  patch.points_are_available=true;
121  // copy points to patch.data
122  for (unsigned int i=0; i<DoFHandlerType::space_dimension; ++i)
123  for (unsigned int q=0; q<n_q_points; ++q)
124  patch.data(patch.data.size(0)-DoFHandlerType::space_dimension+i,q)=q_points[q][i];
125  }
126  else
127  {
128  patch.data.reinit(scratch_data.n_datasets, n_q_points);
129  patch.points_are_available = false;
130  }
131 
132 
133  // counter for data records
134  unsigned int offset=0;
135 
136  // first fill dof_data
137  for (unsigned int dataset=0; dataset<this->dof_data.size(); ++dataset)
138  {
140  = scratch_data.get_present_fe_values (dataset);
141  const unsigned int n_components =
142  this_fe_patch_values.get_fe().n_components();
143 
144  const DataPostprocessor<DoFHandlerType::space_dimension> *postprocessor=this->dof_data[dataset]->postprocessor;
145 
146  if (postprocessor != 0)
147  {
148  // we have to postprocess the data, so determine, which fields
149  // have to be updated
150  const UpdateFlags update_flags=postprocessor->get_needed_update_flags();
151  if (n_components == 1)
152  {
153  // at each point there is only one component of value,
154  // gradient etc.
155  if (update_flags & update_values)
156  this->dof_data[dataset]->get_function_values (this_fe_patch_values,
157  scratch_data.patch_values);
158  if (update_flags & update_gradients)
159  this->dof_data[dataset]->get_function_gradients (this_fe_patch_values,
160  scratch_data.patch_gradients);
161  if (update_flags & update_hessians)
162  this->dof_data[dataset]->get_function_hessians (this_fe_patch_values,
163  scratch_data.patch_hessians);
164 
165  if (update_flags & update_quadrature_points)
166  scratch_data.patch_evaluation_points = this_fe_patch_values.get_quadrature_points();
167 
168 
169  std::vector<Point<DoFHandlerType::space_dimension> > dummy_normals;
170  postprocessor->
171  compute_derived_quantities_scalar(scratch_data.patch_values,
172  scratch_data.patch_gradients,
173  scratch_data.patch_hessians,
174  dummy_normals,
175  scratch_data.patch_evaluation_points,
176  scratch_data.postprocessed_values[dataset]);
177  }
178  else
179  {
180  scratch_data.resize_system_vectors (n_components);
181 
182  // at each point there is a vector valued function and its
183  // derivative...
184  if (update_flags & update_values)
185  this->dof_data[dataset]->get_function_values (this_fe_patch_values,
186  scratch_data.patch_values_system);
187  if (update_flags & update_gradients)
188  this->dof_data[dataset]->get_function_gradients (this_fe_patch_values,
189  scratch_data.patch_gradients_system);
190  if (update_flags & update_hessians)
191  this->dof_data[dataset]->get_function_hessians (this_fe_patch_values,
192  scratch_data.patch_hessians_system);
193 
194  if (update_flags & update_quadrature_points)
195  scratch_data.patch_evaluation_points = this_fe_patch_values.get_quadrature_points();
196 
197  std::vector<Point<DoFHandlerType::space_dimension> > dummy_normals;
198 
199  postprocessor->
200  compute_derived_quantities_vector(scratch_data.patch_values_system,
201  scratch_data.patch_gradients_system,
202  scratch_data.patch_hessians_system,
203  dummy_normals,
204  scratch_data.patch_evaluation_points,
205  scratch_data.postprocessed_values[dataset]);
206  }
207 
208  for (unsigned int q=0; q<n_q_points; ++q)
209  for (unsigned int component=0;
210  component<this->dof_data[dataset]->n_output_variables;
211  ++component)
212  patch.data(offset+component,q)
213  = scratch_data.postprocessed_values[dataset][q](component);
214  }
215  else
216  // now we use the given data vector without modifications. again,
217  // we treat single component functions separately for efficiency
218  // reasons.
219  if (n_components == 1)
220  {
221  this->dof_data[dataset]->get_function_values (this_fe_patch_values,
222  scratch_data.patch_values);
223  for (unsigned int q=0; q<n_q_points; ++q)
224  patch.data(offset,q) = scratch_data.patch_values[q];
225  }
226  else
227  {
228  scratch_data.resize_system_vectors(n_components);
229  this->dof_data[dataset]->get_function_values (this_fe_patch_values,
230  scratch_data.patch_values_system);
231  for (unsigned int component=0; component<n_components;
232  ++component)
233  for (unsigned int q=0; q<n_q_points; ++q)
234  patch.data(offset+component,q) =
235  scratch_data.patch_values_system[q](component);
236  }
237  // increment the counter for the actual data record
238  offset+=this->dof_data[dataset]->n_output_variables;
239  }
240 
241  // then do the cell data. only compute the number of a cell if needed;
242  // also make sure that we only access cell data if the
243  // first_cell/next_cell functions only return active cells
244  if (this->cell_data.size() != 0)
245  {
246  Assert (!cell_and_index->first->has_children(), ExcNotImplemented());
247 
248  for (unsigned int dataset=0; dataset<this->cell_data.size(); ++dataset)
249  {
250  const double value
251  = this->cell_data[dataset]->get_cell_data_value (cell_and_index->second);
252  for (unsigned int q=0; q<n_q_points; ++q)
253  patch.data(offset+dataset,q) = value;
254  }
255  }
256  }
257 
258 
259  for (unsigned int f=0; f<GeometryInfo<DoFHandlerType::dimension>::faces_per_cell; ++f)
260  {
261  // let's look up whether the neighbor behind that face is noted in the
262  // table of cells which we treat. this can only happen if the neighbor
263  // exists, and is on the same level as this cell, but it may also happen
264  // that the neighbor is not a member of the range of cells over which we
265  // loop, in which case the respective entry in the
266  // cell_to_patch_index_map will have the value no_neighbor. (note that
267  // since we allocated only as much space in this array as the maximum
268  // index of the cells we loop over, not every neighbor may have its
269  // space in it, so we have to assume that it is extended by values
270  // no_neighbor)
271  if (cell_and_index->first->at_boundary(f)
272  ||
273  (cell_and_index->first->neighbor(f)->level() != cell_and_index->first->level()))
274  {
276  continue;
277  }
278 
279  const cell_iterator neighbor = cell_and_index->first->neighbor(f);
280  Assert (static_cast<unsigned int>(neighbor->level()) <
281  scratch_data.cell_to_patch_index_map->size(),
282  ExcInternalError());
283  if ((static_cast<unsigned int>(neighbor->index()) >=
284  (*scratch_data.cell_to_patch_index_map)[neighbor->level()].size())
285  ||
286  ((*scratch_data.cell_to_patch_index_map)[neighbor->level()][neighbor->index()]
287  ==
289  {
291  continue;
292  }
293 
294  // now, there is a neighbor, so get its patch number and set it for the
295  // neighbor index
296  patch.neighbors[f]
297  = (*scratch_data.cell_to_patch_index_map)[neighbor->level()][neighbor->index()];
298  }
299 
300  const unsigned int patch_idx =
301  (*scratch_data.cell_to_patch_index_map)[cell_and_index->first->level()][cell_and_index->first->index()];
302  // did we mess up the indices?
303  Assert(patch_idx < patches.size(), ExcInternalError());
304  patch.patch_index = patch_idx;
305 
306  // Put the patch into the patches vector. instead of copying the data,
307  // simply swap the contents to avoid the penalty of writing into another
308  // processor's memory
309  patches[patch_idx].swap (patch);
310 }
311 
312 
313 
314 template <int dim, typename DoFHandlerType>
315 void DataOut<dim,DoFHandlerType>::build_patches (const unsigned int n_subdivisions)
316 {
318  n_subdivisions, no_curved_cells);
319 }
320 
321 
322 
323 template <int dim, typename DoFHandlerType>
326  const unsigned int n_subdivisions_,
327  const CurvedCellRegion curved_region)
328 {
329  // Check consistency of redundant template parameter
330  Assert (dim==DoFHandlerType::dimension, ExcDimensionMismatch(dim, DoFHandlerType::dimension));
331 
332  Assert (this->triangulation != 0,
333  Exceptions::DataOut::ExcNoTriangulationSelected());
334 
335  const unsigned int n_subdivisions = (n_subdivisions_ != 0)
336  ? n_subdivisions_
337  : this->default_subdivisions;
338  Assert (n_subdivisions >= 1,
339  Exceptions::DataOut::ExcInvalidNumberOfSubdivisions(n_subdivisions));
340 
341  // First count the cells we want to create patches of. Also fill the object
342  // that maps the cell indices to the patch numbers, as this will be needed
343  // for generation of neighborship information.
344  // Note, there is a confusing mess of different indices here at play:
345  // patch_index - the index of a patch in all_cells
346  // cell->index - only unique on each level, used in cell_to_patch_index_map
347  // active_index - index for a cell when counting from begin_active() using ++cell
348  // cell_index - unique index of a cell counted using next_locally_owned_cell()
349  // starting from first_locally_owned_cell()
350  //
351  // It turns out that we create one patch for each selected cell, so patch_index==cell_index.
352  //
353  // will be cell_to_patch_index_map[cell->level][cell->index] = patch_index
354  std::vector<std::vector<unsigned int> > cell_to_patch_index_map;
355  cell_to_patch_index_map.resize (this->triangulation->n_levels());
356  for (unsigned int l=0; l<this->triangulation->n_levels(); ++l)
357  {
358  // max_index is the largest cell->index on level l
359  unsigned int max_index = 0;
360  for (cell_iterator cell=first_locally_owned_cell(); cell != this->triangulation->end();
361  cell = next_locally_owned_cell(cell))
362  if (static_cast<unsigned int>(cell->level()) == l)
363  max_index = std::max (max_index,
364  static_cast<unsigned int>(cell->index()));
365 
366  cell_to_patch_index_map[l].resize (max_index+1,
367  ::DataOutBase::Patch<DoFHandlerType::dimension,
368  DoFHandlerType::space_dimension>::no_neighbor);
369  }
370 
371  // will be all_cells[patch_index] = pair(cell, active_index)
372  std::vector<std::pair<cell_iterator, unsigned int> > all_cells;
373  {
374  // important: we need to compute the active_index of the cell in the range
375  // 0..n_active_cells() because this is where we need to look up cell
376  // data from (cell data vectors do not have the length distance computed by
377  // first_locally_owned_cell/next_locally_owned_cell because this might skip
378  // some values (FilteredIterator).
379  active_cell_iterator active_cell = this->triangulation->begin_active();
380  unsigned int active_index = 0;
381  cell_iterator cell = first_locally_owned_cell();
382  for (; cell != this->triangulation->end();
383  cell = next_locally_owned_cell(cell))
384  {
385  // move forward until active_cell points at the cell (cell) we are looking
386  // at to compute the current active_index
387  while (active_cell!=this->triangulation->end()
388  && cell->active()
389  && active_cell_iterator(cell) != active_cell)
390  {
391  ++active_cell;
392  ++active_index;
393  }
394 
395  Assert (static_cast<unsigned int>(cell->level()) <
396  cell_to_patch_index_map.size(),
397  ExcInternalError());
398  Assert (static_cast<unsigned int>(cell->index()) <
399  cell_to_patch_index_map[cell->level()].size(),
400  ExcInternalError());
401  Assert (active_index < this->triangulation->n_active_cells(),
402  ExcInternalError());
403  cell_to_patch_index_map[cell->level()][cell->index()] = all_cells.size();
404 
405  all_cells.push_back (std::make_pair(cell, active_index));
406  }
407  }
408 
409  this->patches.clear ();
410  this->patches.resize(all_cells.size());
411 
412  // now create a default object for the WorkStream object to work with
413  unsigned int n_datasets=this->cell_data.size();
414  for (unsigned int i=0; i<this->dof_data.size(); ++i)
415  n_datasets += this->dof_data[i]->n_output_variables;
416 
417  std::vector<unsigned int> n_postprocessor_outputs (this->dof_data.size());
418  for (unsigned int dataset=0; dataset<this->dof_data.size(); ++dataset)
419  if (this->dof_data[dataset]->postprocessor)
420  n_postprocessor_outputs[dataset] = this->dof_data[dataset]->n_output_variables;
421  else
422  n_postprocessor_outputs[dataset] = 0;
423 
424  const CurvedCellRegion curved_cell_region
425  = (n_subdivisions<2 ? no_curved_cells : curved_region);
426 
427  UpdateFlags update_flags = update_values;
428  if (curved_cell_region != no_curved_cells)
429  update_flags |= update_quadrature_points;
430 
431  for (unsigned int i=0; i<this->dof_data.size(); ++i)
432  if (this->dof_data[i]->postprocessor)
433  update_flags |= this->dof_data[i]->postprocessor->get_needed_update_flags();
434  // perhaps update_normal_vectors is present, which would only be useful on
435  // faces, but we may not use it here.
436  Assert (!(update_flags & update_normal_vectors),
437  ExcMessage("The update of normal vectors may not be requested for evaluation of "
438  "data on cells via DataPostprocessor."));
439 
441  thread_data (n_datasets, n_subdivisions,
442  n_postprocessor_outputs,
443  mapping,
444  this->get_finite_elements(),
445  update_flags,
446  cell_to_patch_index_map);
447 
448  // now build the patches in parallel
449  if (all_cells.size() > 0)
450  WorkStream::run (&all_cells[0],
451  &all_cells[0]+all_cells.size(),
453  this,
454  std_cxx11::_1,
455  std_cxx11::_2,
456  /* no std_cxx11::_3, since this function doesn't actually need a
457  copy data object -- it just writes everything right into the
458  output array */
459  n_subdivisions,
460  curved_cell_region,
461  std_cxx11::ref(this->patches)),
462  // no copy-local-to-global function needed here
463  std_cxx11::function<void (const int &)>(),
464  thread_data,
465  /* dummy CopyData object = */ 0,
466  // experimenting shows that we can make things run a bit
467  // faster if we increase the number of cells we work on
468  // per item (i.e., WorkStream's chunk_size argument,
469  // about 10% improvement) and the items in flight at any
470  // given time (another 5% on the testcase discussed in
471  // @ref workstream_paper, on 32 cores) and if
473  64);
474 }
475 
476 
477 
478 template <int dim, typename DoFHandlerType>
481 {
482  return this->triangulation->begin_active ();
483 }
484 
485 
486 
487 template <int dim, typename DoFHandlerType>
491 {
492  // convert the iterator to an active_iterator and advance this to the next
493  // active cell
495  active_cell_iterator active_cell = cell;
496  ++active_cell;
497  return active_cell;
498 }
499 
500 
501 
502 template <int dim, typename DoFHandlerType>
505 {
507  cell = first_cell();
508 
509  // skip cells if the current one has no children (is active) and is a ghost
510  // or artificial cell
511  while ((cell != this->triangulation->end()) &&
512  (cell->has_children() == false) &&
513  !cell->is_locally_owned())
514  cell = next_cell(cell);
515 
516  return cell;
517 }
518 
519 
520 
521 template <int dim, typename DoFHandlerType>
525 {
527  cell = next_cell(old_cell);
528  while ((cell != this->triangulation->end()) &&
529  (cell->has_children() == false) &&
530  !cell->is_locally_owned())
531  cell = next_cell(cell);
532  return cell;
533 }
534 
535 
536 // explicit instantiations
537 #include "data_out.inst"
538 
539 DEAL_II_NAMESPACE_CLOSE
Shape function values.
static const unsigned int invalid_unsigned_int
Definition: types.h:164
void build_one_patch(const std::pair< cell_iterator, unsigned int > *cell_and_index, internal::DataOut::ParallelData< DoFHandlerType::dimension, DoFHandlerType::space_dimension > &scratch_data, const unsigned int n_subdivisions, const CurvedCellRegion curved_cell_region, std::vector< DataOutBase::Patch< DoFHandlerType::dimension, DoFHandlerType::space_dimension > > &patches)
Definition: data_out.cc:66
::ExceptionBase & ExcMessage(std::string arg1)
virtual cell_iterator next_cell(const cell_iterator &cell)
Definition: data_out.cc:490
Transformed quadrature points.
Point< spacedim > vertices[GeometryInfo< dim >::vertices_per_cell]
virtual void build_patches(const unsigned int n_subdivisions=0)
Definition: data_out.cc:315
unsigned int patch_index
const std::vector< Point< spacedim > > & get_quadrature_points() const
virtual cell_iterator first_locally_owned_cell()
Definition: data_out.cc:504
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
#define Assert(cond, exc)
Definition: exceptions.h:294
UpdateFlags
Abstract base class for mapping classes.
Definition: dof_tools.h:52
virtual cell_iterator first_cell()
Definition: data_out.cc:480
unsigned int n_subdivisions
DataOut_DoFData< DoFHandlerType, DoFHandlerType::dimension, DoFHandlerType::space_dimension >::cell_iterator cell_iterator
Definition: data_out.h:166
Second derivatives of shape functions.
const unsigned int n_quadrature_points
Definition: fe_values.h:1457
static const unsigned int space_dim
Shape function gradients.
Normal vectors.
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:1119
const FiniteElement< dim, spacedim > & get_fe() const
static unsigned int n_threads()
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
unsigned int size(const unsigned int i) const
CurvedCellRegion
Definition: data_out.h:174
virtual cell_iterator next_locally_owned_cell(const cell_iterator &cell)
Definition: data_out.cc:524
unsigned int neighbors[dim > 0?GeometryInfo< dim >::faces_per_cell:1]
Table< 2, float > data
virtual UpdateFlags get_needed_update_flags() const =0