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