Reference documentation for deal.II version GIT c415589cf0 2022-08-14 18:50:02+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_stack.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2021 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 VectorType vector_type)
86 {
87  std::vector<std::string> names;
88  names.push_back(name);
89  declare_data_vector(names, vector_type);
90 }
91 
92 
93 template <int dim, int spacedim>
94 void
96  const std::vector<std::string> &names,
97  const VectorType vector_type)
98 {
99 #ifdef DEBUG
100  // make sure this function is
101  // not called after some parameter
102  // values have already been
103  // processed
104  Assert(patches.size() == 0, ExcDataAlreadyAdded());
105 
106  // also make sure that no name is
107  // used twice
108  for (const auto &name : names)
109  {
110  for (const auto &data_set : dof_data)
111  for (const auto &data_set_name : data_set.names)
112  Assert(name != data_set_name, ExcNameAlreadyUsed(name));
113 
114  for (const auto &data_set : cell_data)
115  for (const auto &data_set_name : data_set.names)
116  Assert(name != data_set_name, ExcNameAlreadyUsed(name));
117  }
118 #endif
119 
120  switch (vector_type)
121  {
122  case dof_vector:
123  dof_data.emplace_back();
124  dof_data.back().names = names;
125  break;
126 
127  case cell_vector:
128  cell_data.emplace_back();
129  cell_data.back().names = names;
130  break;
131  }
132 }
133 
134 
135 template <int dim, int spacedim>
136 template <typename number>
137 void
139  const std::string & name)
140 {
141  const unsigned int n_components = dof_handler->get_fe(0).n_components();
142 
143  std::vector<std::string> names;
144  // if only one component or vector
145  // is cell vector: we only need one
146  // name
147  if ((n_components == 1) ||
148  (vec.size() == dof_handler->get_triangulation().n_active_cells()))
149  {
150  names.resize(1, name);
151  }
152  else
153  // otherwise append _i to the
154  // given name
155  {
156  names.resize(n_components);
157  for (unsigned int i = 0; i < n_components; ++i)
158  {
159  std::ostringstream namebuf;
160  namebuf << '_' << i;
161  names[i] = name + namebuf.str();
162  }
163  }
164 
165  add_data_vector(vec, names);
166 }
167 
168 
169 template <int dim, int spacedim>
170 template <typename number>
171 void
173  const Vector<number> & vec,
174  const std::vector<std::string> &names)
175 {
176  Assert(dof_handler != nullptr,
178  // either cell data and one name,
179  // or dof data and n_components names
180  Assert(((vec.size() == dof_handler->get_triangulation().n_active_cells()) &&
181  (names.size() == 1)) ||
182  ((vec.size() == dof_handler->n_dofs()) &&
183  (names.size() == dof_handler->get_fe(0).n_components())),
185  names.size(), dof_handler->get_fe(0).n_components()));
186  for (const auto &name : names)
187  {
188  (void)name;
189  Assert(name.find_first_not_of("abcdefghijklmnopqrstuvwxyz"
190  "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
191  "0123456789_<>()") == std::string::npos,
193  name,
194  name.find_first_not_of("abcdefghijklmnopqrstuvwxyz"
195  "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
196  "0123456789_<>()")));
197  }
198 
199  if (vec.size() == dof_handler->n_dofs())
200  {
201  typename std::vector<DataVector>::iterator data_vector = dof_data.begin();
202  for (; data_vector != dof_data.end(); ++data_vector)
203  if (data_vector->names == names)
204  {
205  data_vector->data.reinit(vec.size());
206  std::copy(vec.begin(), vec.end(), data_vector->data.begin());
207  return;
208  }
209 
210  // ok. not found. there is a
211  // slight chance that
212  // n_dofs==n_cells, so only
213  // bomb out if the next if
214  // statement will not be run
215  if (dof_handler->n_dofs() !=
216  dof_handler->get_triangulation().n_active_cells())
217  Assert(false, ExcVectorNotDeclared(names[0]));
218  }
219 
220  // search cell data
221  if ((vec.size() != dof_handler->n_dofs()) ||
222  (dof_handler->n_dofs() ==
223  dof_handler->get_triangulation().n_active_cells()))
224  {
225  typename std::vector<DataVector>::iterator data_vector =
226  cell_data.begin();
227  for (; data_vector != cell_data.end(); ++data_vector)
228  if (data_vector->names == names)
229  {
230  data_vector->data.reinit(vec.size());
231  std::copy(vec.begin(), vec.end(), data_vector->data.begin());
232  return;
233  }
234  Assert(false, ExcVectorNotDeclared(names[0]));
235  }
236 
237  // we have either return or Assert
238  // statements above, so shouldn't
239  // get here!
240  Assert(false, ExcInternalError());
241 }
242 
243 
244 template <int dim, int spacedim>
245 void
246 DataOutStack<dim, spacedim>::build_patches(const unsigned int nnnn_subdivisions)
247 {
248  // this is mostly copied from the
249  // DataOut class
250  unsigned int n_subdivisions =
251  (nnnn_subdivisions != 0) ? nnnn_subdivisions : this->default_subdivisions;
252 
253  Assert(n_subdivisions >= 1,
255  n_subdivisions));
256  Assert(dof_handler != nullptr,
258 
259  this->validate_dataset_names();
260 
261  const unsigned int n_components = dof_handler->get_fe(0).n_components();
262  const unsigned int n_datasets =
263  dof_data.size() * n_components + cell_data.size();
264 
265  // first count the cells we want to
266  // create patches of and make sure
267  // there is enough memory for that
268  unsigned int n_patches = 0;
270  dof_handler->begin_active();
271  cell != dof_handler->end();
272  ++cell)
273  ++n_patches;
274 
275 
276  // before we start the loop:
277  // create a quadrature rule that
278  // actually has the points on this
279  // patch, and an object that
280  // extracts the data on each
281  // cell to these points
282  QTrapezoid<1> q_trapez;
283  QIterated<dim> patch_points(q_trapez, n_subdivisions);
284 
285  // create collection objects from
286  // single quadratures,
287  // and finite elements. if we have
288  // an hp-DoFHandler,
289  // dof_handler.get_fe() returns a
290  // collection of which we do a
291  // shallow copy instead
292  const hp::QCollection<dim> q_collection(patch_points);
293  const hp::FECollection<dim> &fe_collection = dof_handler->get_fe_collection();
294 
295  hp::FEValues<dim> x_fe_patch_values(fe_collection,
296  q_collection,
297  update_values);
298 
299  const unsigned int n_q_points = patch_points.size();
300  std::vector<double> patch_values(n_q_points);
301  std::vector<Vector<double>> patch_values_system(n_q_points,
302  Vector<double>(n_components));
303 
304  // add the required number of
305  // patches. first initialize a template
306  // patch with n_q_points (in the plane
307  // of the cells) times n_subdivisions+1 (in
308  // the time direction) points
310  default_patch.n_subdivisions = n_subdivisions;
311  default_patch.reference_cell = ReferenceCells::get_hypercube<dim + 1>();
312  default_patch.data.reinit(n_datasets, n_q_points * (n_subdivisions + 1));
313  patches.insert(patches.end(), n_patches, default_patch);
314 
315  // now loop over all cells and
316  // actually create the patches
317  typename std::vector<
319  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),
354  patch->vertices[1] = Point<dim + 1>(cell->vertex(1)(0),
355  cell->vertex(1)(1),
357  patch->vertices[2] = Point<dim + 1>(cell->vertex(2)(0),
358  cell->vertex(2)(1),
360  patch->vertices[3] = Point<dim + 1>(cell->vertex(3)(0),
361  cell->vertex(3)(1),
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 {
461 }
462 
463 
464 
465 template <int dim, int spacedim>
466 const std::vector<
470 {
471  return patches;
472 }
473 
474 
475 
476 template <int dim, int spacedim>
477 std::vector<std::string>
479 {
480  std::vector<std::string> names;
481  for (typename std::vector<DataVector>::const_iterator dataset =
482  dof_data.begin();
483  dataset != dof_data.end();
484  ++dataset)
485  names.insert(names.end(), dataset->names.begin(), dataset->names.end());
486  for (typename std::vector<DataVector>::const_iterator dataset =
487  cell_data.begin();
488  dataset != cell_data.end();
489  ++dataset)
490  names.insert(names.end(), dataset->names.begin(), dataset->names.end());
491 
492  return names;
493 }
494 
495 
496 
497 // explicit instantiations
498 #include "data_out_stack.inst"
499 
500 
std::size_t memory_consumption() const
void declare_data_vector(const std::string &name, const VectorType vector_type)
virtual std::vector< std::string > get_dataset_names() const override
std::vector< DataVector > cell_data
void new_parameter_value(const double parameter_value, const double parameter_step)
std::vector< DataVector > dof_data
SmartPointer< const DoFHandler< dim, spacedim >, DataOutStack< dim, spacedim > > dof_handler
void finish_parameter_value()
void attach_dof_handler(const DoFHandler< dim, spacedim > &dof_handler)
double parameter_step
void build_patches(const unsigned int n_subdivisions=0)
std::vector<::DataOutBase::Patch< patch_dim, patch_spacedim > > patches
virtual const std::vector< ::DataOutBase::Patch< DataOutStack< dim, spacedim >::patch_dim, DataOutStack< dim, spacedim >::patch_spacedim > > & get_patches() const override
void add_data_vector(const Vector< number > &vec, const std::string &name)
void get_function_values(const InputVector &fe_function, std::vector< typename InputVector::value_type > &values) const
Definition: fe_values.cc:3356
Definition: point.h:111
unsigned int size() const
Definition: vector.h:110
iterator end()
size_type size() const
iterator begin()
const FEValuesType & get_present_fe_values() const
Definition: fe_values.h:693
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:294
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNoDoFHandlerSelected()
static ::ExceptionBase & ExcInvalidNumberOfSubdivisions(int arg1)
static ::ExceptionBase & ExcInvalidCharacter(std::string arg1, size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcNameAlreadyUsed(std::string arg1)
static ::ExceptionBase & ExcVectorNotDeclared(std::string arg1)
static ::ExceptionBase & ExcInvalidNumberOfNames(int arg1, int arg2)
static ::ExceptionBase & ExcDataAlreadyAdded()
static ::ExceptionBase & ExcDataNotCleared()
@ update_values
Shape function values.
std::enable_if_t< std::is_fundamental< T >::value, std::size_t > memory_consumption(const T &t)
void copy(const T *begin, const T *end, U *dest)
ReferenceCell reference_cell
Table< 2, float > data
unsigned int n_subdivisions
std::size_t memory_consumption() const
std::vector< std::string > names