Reference documentation for deal.II version Git 7026f387cc 2019-10-15 14:19:01 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
data_out_stack.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 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/memory_consumption.h>
17 #include <deal.II/base/quadrature_lib.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_iterator.h>
27 
28 #include <deal.II/hp/fe_values.h>
29 
30 #include <deal.II/lac/block_vector.h>
31 #include <deal.II/lac/vector.h>
32 
33 #include <deal.II/numerics/data_out_stack.h>
34 
35 #include <sstream>
36 
37 DEAL_II_NAMESPACE_OPEN
38 
39 
40 template <int dim, int spacedim, typename DoFHandlerType>
41 std::size_t
43  const
44 {
47 }
48 
49 
50 
51 template <int dim, int spacedim, typename DoFHandlerType>
52 void
54  const double p,
55  const double dp)
56 {
57  parameter = p;
58  parameter_step = dp;
59 
60  // check whether the user called finish_parameter_value() at the end of the
61  // previous parameter step
62  //
63  // this is to prevent serious waste of memory
64  for (typename std::vector<DataVector>::const_iterator i = dof_data.begin();
65  i != dof_data.end();
66  ++i)
67  Assert(i->data.size() == 0, ExcDataNotCleared());
68  for (typename std::vector<DataVector>::const_iterator i = cell_data.begin();
69  i != cell_data.end();
70  ++i)
71  Assert(i->data.size() == 0, ExcDataNotCleared());
72 }
73 
74 
75 template <int dim, int spacedim, typename DoFHandlerType>
76 void
78  const DoFHandlerType &dof)
79 {
80  // Check consistency of redundant
81  // template parameter
82  Assert(dim == DoFHandlerType::dimension,
83  ExcDimensionMismatch(dim, DoFHandlerType::dimension));
84 
85  dof_handler = &dof;
86 }
87 
88 
89 template <int dim, int spacedim, typename DoFHandlerType>
90 void
92  const std::string &name,
93  const VectorType vector_type)
94 {
95  std::vector<std::string> names;
96  names.push_back(name);
97  declare_data_vector(names, vector_type);
98 }
99 
100 
101 template <int dim, int spacedim, typename DoFHandlerType>
102 void
104  const std::vector<std::string> &names,
105  const VectorType vector_type)
106 {
107 #ifdef DEBUG
108  // make sure this function is
109  // not called after some parameter
110  // values have already been
111  // processed
112  Assert(patches.size() == 0, ExcDataAlreadyAdded());
113 
114  // also make sure that no name is
115  // used twice
116  for (const auto &name : names)
117  {
118  for (const auto &data_set : dof_data)
119  for (const auto &data_set_name : data_set.names)
120  Assert(name != data_set_name, ExcNameAlreadyUsed(name));
121 
122  for (const auto &data_set : cell_data)
123  for (const auto &data_set_name : data_set.names)
124  Assert(name != data_set_name, ExcNameAlreadyUsed(name));
125  }
126 #endif
127 
128  switch (vector_type)
129  {
130  case dof_vector:
131  dof_data.emplace_back();
132  dof_data.back().names = names;
133  break;
134 
135  case cell_vector:
136  cell_data.emplace_back();
137  cell_data.back().names = names;
138  break;
139  }
140 }
141 
142 
143 template <int dim, int spacedim, typename DoFHandlerType>
144 template <typename number>
145 void
147  const Vector<number> &vec,
148  const std::string & name)
149 {
150  const unsigned int n_components = dof_handler->get_fe(0).n_components();
151 
152  std::vector<std::string> names;
153  // if only one component or vector
154  // is cell vector: we only need one
155  // name
156  if ((n_components == 1) ||
157  (vec.size() == dof_handler->get_triangulation().n_active_cells()))
158  {
159  names.resize(1, name);
160  }
161  else
162  // otherwise append _i to the
163  // given name
164  {
165  names.resize(n_components);
166  for (unsigned int i = 0; i < n_components; ++i)
167  {
168  std::ostringstream namebuf;
169  namebuf << '_' << i;
170  names[i] = name + namebuf.str();
171  }
172  }
173 
174  add_data_vector(vec, names);
175 }
176 
177 
178 template <int dim, int spacedim, typename DoFHandlerType>
179 template <typename number>
180 void
182  const Vector<number> & vec,
183  const std::vector<std::string> &names)
184 {
185  Assert(dof_handler != nullptr,
187  // either cell data and one name,
188  // or dof data and n_components names
189  Assert(((vec.size() == dof_handler->get_triangulation().n_active_cells()) &&
190  (names.size() == 1)) ||
191  ((vec.size() == dof_handler->n_dofs()) &&
192  (names.size() == dof_handler->get_fe(0).n_components())),
194  names.size(), dof_handler->get_fe(0).n_components()));
195  for (const auto &name : names)
196  {
197  (void)name;
198  Assert(name.find_first_not_of("abcdefghijklmnopqrstuvwxyz"
199  "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
200  "0123456789_<>()") == std::string::npos,
202  name,
203  name.find_first_not_of("abcdefghijklmnopqrstuvwxyz"
204  "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
205  "0123456789_<>()")));
206  }
207 
208  if (vec.size() == dof_handler->n_dofs())
209  {
210  typename std::vector<DataVector>::iterator data_vector = dof_data.begin();
211  for (; data_vector != dof_data.end(); ++data_vector)
212  if (data_vector->names == names)
213  {
214  data_vector->data.reinit(vec.size());
215  std::copy(vec.begin(), vec.end(), data_vector->data.begin());
216  return;
217  }
218 
219  // ok. not found. there is a
220  // slight chance that
221  // n_dofs==n_cells, so only
222  // bomb out if the next if
223  // statement will not be run
224  if (dof_handler->n_dofs() !=
225  dof_handler->get_triangulation().n_active_cells())
226  Assert(false, ExcVectorNotDeclared(names[0]));
227  }
228 
229  // search cell data
230  if ((vec.size() != dof_handler->n_dofs()) ||
231  (dof_handler->n_dofs() ==
232  dof_handler->get_triangulation().n_active_cells()))
233  {
234  typename std::vector<DataVector>::iterator data_vector =
235  cell_data.begin();
236  for (; data_vector != cell_data.end(); ++data_vector)
237  if (data_vector->names == names)
238  {
239  data_vector->data.reinit(vec.size());
240  std::copy(vec.begin(), vec.end(), data_vector->data.begin());
241  return;
242  }
243  Assert(false, ExcVectorNotDeclared(names[0]));
244  }
245 
246  // we have either return or Assert
247  // statements above, so shouldn't
248  // get here!
249  Assert(false, ExcInternalError());
250 }
251 
252 
253 template <int dim, int spacedim, typename DoFHandlerType>
254 void
256  const unsigned int nnnn_subdivisions)
257 {
258  // this is mostly copied from the
259  // DataOut class
260  unsigned int n_subdivisions =
261  (nnnn_subdivisions != 0) ? nnnn_subdivisions : this->default_subdivisions;
262 
263  Assert(n_subdivisions >= 1,
265  n_subdivisions));
266  Assert(dof_handler != nullptr,
268 
269  this->validate_dataset_names();
270 
271  const unsigned int n_components = dof_handler->get_fe(0).n_components();
272  const unsigned int n_datasets =
273  dof_data.size() * n_components + cell_data.size();
274 
275  // first count the cells we want to
276  // create patches of and make sure
277  // there is enough memory for that
278  unsigned int n_patches = 0;
279  for (typename DoFHandlerType::active_cell_iterator cell =
280  dof_handler->begin_active();
281  cell != dof_handler->end();
282  ++cell)
283  ++n_patches;
284 
285 
286  // before we start the loop:
287  // create a quadrature rule that
288  // actually has the points on this
289  // patch, and an object that
290  // extracts the data on each
291  // cell to these points
292  QTrapez<1> q_trapez;
293  QIterated<dim> patch_points(q_trapez, n_subdivisions);
294 
295  // create collection objects from
296  // single quadratures,
297  // and finite elements. if we have
298  // an hp DoFHandler,
299  // dof_handler.get_fe() returns a
300  // collection of which we do a
301  // shallow copy instead
302  const hp::QCollection<dim> q_collection(patch_points);
303  const hp::FECollection<dim> &fe_collection = dof_handler->get_fe_collection();
304 
305  hp::FEValues<dim> x_fe_patch_values(fe_collection,
306  q_collection,
307  update_values);
308 
309  const unsigned int n_q_points = patch_points.size();
310  std::vector<double> patch_values(n_q_points);
311  std::vector<Vector<double>> patch_values_system(n_q_points,
312  Vector<double>(n_components));
313 
314  // add the required number of
315  // patches. first initialize a template
316  // patch with n_q_points (in the plane
317  // of the cells) times n_subdivisions+1 (in
318  // the time direction) points
320  default_patch.n_subdivisions = n_subdivisions;
321  default_patch.data.reinit(n_datasets, n_q_points * (n_subdivisions + 1));
322  patches.insert(patches.end(), n_patches, default_patch);
323 
324  // now loop over all cells and
325  // actually create the patches
326  typename std::vector<::DataOutBase::Patch<dim + 1, dim + 1>>::iterator
327  patch = patches.begin() + (patches.size() - n_patches);
328  unsigned int cell_number = 0;
329  for (typename DoFHandlerType::active_cell_iterator cell =
330  dof_handler->begin_active();
331  cell != dof_handler->end();
332  ++cell, ++patch, ++cell_number)
333  {
334  Assert(cell->is_locally_owned(), ExcNotImplemented());
335 
336  Assert(patch != patches.end(), ExcInternalError());
337 
338  // first fill in the vertices of the patch
339 
340  // Patches are organized such
341  // that the parameter direction
342  // is the last
343  // coordinate. Thus, vertices
344  // are two copies of the space
345  // patch, one at parameter-step
346  // and one at parameter.
347  switch (dim)
348  {
349  case 1:
350  patch->vertices[0] =
351  Point<dim + 1>(cell->vertex(0)(0), parameter - parameter_step);
352  patch->vertices[1] =
353  Point<dim + 1>(cell->vertex(1)(0), parameter - parameter_step);
354  patch->vertices[2] = Point<dim + 1>(cell->vertex(0)(0), parameter);
355  patch->vertices[3] = Point<dim + 1>(cell->vertex(1)(0), parameter);
356  break;
357 
358  case 2:
359  patch->vertices[0] = Point<dim + 1>(cell->vertex(0)(0),
360  cell->vertex(0)(1),
362  patch->vertices[1] = Point<dim + 1>(cell->vertex(1)(0),
363  cell->vertex(1)(1),
365  patch->vertices[2] = Point<dim + 1>(cell->vertex(2)(0),
366  cell->vertex(2)(1),
368  patch->vertices[3] = Point<dim + 1>(cell->vertex(3)(0),
369  cell->vertex(3)(1),
371  patch->vertices[4] =
372  Point<dim + 1>(cell->vertex(0)(0), cell->vertex(0)(1), parameter);
373  patch->vertices[5] =
374  Point<dim + 1>(cell->vertex(1)(0), cell->vertex(1)(1), parameter);
375  patch->vertices[6] =
376  Point<dim + 1>(cell->vertex(2)(0), cell->vertex(2)(1), parameter);
377  patch->vertices[7] =
378  Point<dim + 1>(cell->vertex(3)(0), cell->vertex(3)(1), parameter);
379  break;
380 
381  default:
382  Assert(false, ExcNotImplemented());
383  }
384 
385 
386  // now fill in the data values.
387  // note that the required order is
388  // with highest coordinate running
389  // fastest, we need to enter each
390  // value (n_subdivisions+1) times
391  // in succession
392  if (n_datasets > 0)
393  {
394  x_fe_patch_values.reinit(cell);
395  const FEValues<dim> &fe_patch_values =
396  x_fe_patch_values.get_present_fe_values();
397 
398  // first fill dof_data
399  for (unsigned int dataset = 0; dataset < dof_data.size(); ++dataset)
400  {
401  if (n_components == 1)
402  {
403  fe_patch_values.get_function_values(dof_data[dataset].data,
404  patch_values);
405  for (unsigned int i = 0; i < n_subdivisions + 1; ++i)
406  for (unsigned int q = 0; q < n_q_points; ++q)
407  patch->data(dataset, q + n_q_points * i) =
408  patch_values[q];
409  }
410  else
411  // system of components
412  {
413  fe_patch_values.get_function_values(dof_data[dataset].data,
414  patch_values_system);
415  for (unsigned int component = 0; component < n_components;
416  ++component)
417  for (unsigned int i = 0; i < n_subdivisions + 1; ++i)
418  for (unsigned int q = 0; q < n_q_points; ++q)
419  patch->data(dataset * n_components + component,
420  q + n_q_points * i) =
421  patch_values_system[q](component);
422  }
423  }
424 
425  // then do the cell data
426  for (unsigned int dataset = 0; dataset < cell_data.size(); ++dataset)
427  {
428  const double value = cell_data[dataset].data(cell_number);
429  for (unsigned int q = 0; q < n_q_points; ++q)
430  for (unsigned int i = 0; i < n_subdivisions + 1; ++i)
431  patch->data(dataset + dof_data.size() * n_components,
432  q * (n_subdivisions + 1) + i) = value;
433  }
434  }
435  }
436 }
437 
438 
439 template <int dim, int spacedim, typename DoFHandlerType>
440 void
442 {
443  // release lock on dof handler
444  dof_handler = nullptr;
445  for (typename std::vector<DataVector>::iterator i = dof_data.begin();
446  i != dof_data.end();
447  ++i)
448  i->data.reinit(0);
449 
450  for (typename std::vector<DataVector>::iterator i = cell_data.begin();
451  i != cell_data.end();
452  ++i)
453  i->data.reinit(0);
454 }
455 
456 
457 
458 template <int dim, int spacedim, typename DoFHandlerType>
459 std::size_t
461 {
469 }
470 
471 
472 
473 template <int dim, int spacedim, typename DoFHandlerType>
474 const std::vector<::DataOutBase::Patch<dim + 1, dim + 1>> &
476 {
477  return patches;
478 }
479 
480 
481 
482 template <int dim, int spacedim, typename DoFHandlerType>
483 std::vector<std::string>
485 {
486  std::vector<std::string> names;
487  for (typename std::vector<DataVector>::const_iterator dataset =
488  dof_data.begin();
489  dataset != dof_data.end();
490  ++dataset)
491  names.insert(names.end(), dataset->names.begin(), dataset->names.end());
492  for (typename std::vector<DataVector>::const_iterator dataset =
493  cell_data.begin();
494  dataset != cell_data.end();
495  ++dataset)
496  names.insert(names.end(), dataset->names.begin(), dataset->names.end());
497 
498  return names;
499 }
500 
501 
502 
503 // explicit instantiations
504 #include "data_out_stack.inst"
505 
506 
507 DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcNameAlreadyUsed(std::string arg1)
Shape function values.
std::size_t memory_consumption() const
std::vector< DataVector > cell_data
static ::ExceptionBase & ExcDataAlreadyAdded()
void attach_dof_handler(const DoFHandlerType &dof_handler)
void get_function_values(const InputVector &fe_function, std::vector< typename InputVector::value_type > &values) const
Definition: fe_values.cc:3579
std::vector< DataVector > dof_data
unsigned int default_subdivisions
static ::ExceptionBase & ExcInvalidNumberOfSubdivisions(int arg1)
iterator end()
Definition: point.h:110
virtual std::vector< std::string > get_dataset_names() const override
static ::ExceptionBase & ExcDataNotCleared()
SmartPointer< const DoFHandlerType, DataOutStack< dim, spacedim, DoFHandlerType > > dof_handler
std::size_t memory_consumption() const
double parameter_step
std::vector< std::string > names
static ::ExceptionBase & ExcNoDoFHandlerSelected()
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
void declare_data_vector(const std::string &name, const VectorType vector_type)
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcInvalidNumberOfNames(int arg1, int arg2)
virtual const std::vector<::DataOutBase::Patch< dim+1, dim+1 > > & get_patches() const override
static ::ExceptionBase & ExcInvalidCharacter(std::string arg1, size_t arg2)
unsigned int n_subdivisions
void reinit(const TriaIterator< DoFCellAccessor< DoFHandlerType, lda >> cell, const unsigned int q_index=numbers::invalid_unsigned_int, const unsigned int mapping_index=numbers::invalid_unsigned_int, const unsigned int fe_index=numbers::invalid_unsigned_int)
Definition: fe_values.cc:135
iterator begin()
static ::ExceptionBase & ExcVectorNotDeclared(std::string arg1)
unsigned int size() const
void build_patches(const unsigned int n_subdivisions=0)
size_type size() const
std::vector<::DataOutBase::Patch< dim+1, dim+1 > > patches
void validate_dataset_names() const
static ::ExceptionBase & ExcNotImplemented()
void finish_parameter_value()
Table< 2, float > data
void new_parameter_value(const double parameter_value, const double parameter_step)
void add_data_vector(const Vector< number > &vec, const std::string &name)
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
static ::ExceptionBase & ExcInternalError()
const ::FEValues< dim, spacedim > & get_present_fe_values() const
Definition: fe_values.h:569