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