Reference documentation for deal.II version GIT relicensing-233-g802318d791 2024-03-28 20:20: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\}}\)
Loading...
Searching...
No Matches
data_out_stack.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 1999 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
17
20
21#include <deal.II/fe/fe.h>
23
25
27
29#include <deal.II/lac/vector.h>
30
32
33#include <sstream>
34
36
37
38template <int dim, int spacedim>
39std::size_t
45
46
47
48template <int dim, int spacedim>
49void
51 const double dp)
52{
53 parameter = p;
54 parameter_step = dp;
55
56 // check whether the user called finish_parameter_value() at the end of the
57 // previous parameter step
58 //
59 // this is to prevent serious waste of memory
60 for (typename std::vector<DataVector>::const_iterator i = dof_data.begin();
61 i != dof_data.end();
62 ++i)
63 Assert(i->data.size() == 0, ExcDataNotCleared());
64 for (typename std::vector<DataVector>::const_iterator i = cell_data.begin();
65 i != cell_data.end();
66 ++i)
67 Assert(i->data.size() == 0, ExcDataNotCleared());
68}
69
70
71template <int dim, int spacedim>
72void
78
79
80template <int dim, int spacedim>
81void
83 const VectorType vector_type)
84{
85 std::vector<std::string> names;
86 names.push_back(name);
87 declare_data_vector(names, vector_type);
88}
89
90
91template <int dim, int spacedim>
92void
94 const std::vector<std::string> &names,
95 const VectorType vector_type)
96{
97#ifdef DEBUG
98 // make sure this function is
99 // not called after some parameter
100 // values have already been
101 // processed
103
104 // also make sure that no name is
105 // used twice
106 for (const auto &name : names)
107 {
108 for (const auto &data_set : dof_data)
109 for (const auto &data_set_name : data_set.names)
110 Assert(name != data_set_name, ExcNameAlreadyUsed(name));
111
112 for (const auto &data_set : cell_data)
113 for (const auto &data_set_name : data_set.names)
114 Assert(name != data_set_name, ExcNameAlreadyUsed(name));
115 }
116#endif
117
118 switch (vector_type)
119 {
120 case dof_vector:
121 dof_data.emplace_back();
122 dof_data.back().names = names;
123 break;
124
125 case cell_vector:
126 cell_data.emplace_back();
127 cell_data.back().names = names;
128 break;
129 }
130}
131
132
133template <int dim, int spacedim>
134template <typename number>
135void
137 const std::string &name)
138{
139 const unsigned int n_components = dof_handler->get_fe(0).n_components();
140
141 std::vector<std::string> names;
142 // if only one component or vector
143 // is cell vector: we only need one
144 // name
145 if ((n_components == 1) ||
146 (vec.size() == dof_handler->get_triangulation().n_active_cells()))
147 {
148 names.resize(1, name);
149 }
150 else
151 // otherwise append _i to the
152 // given name
153 {
154 names.resize(n_components);
155 for (unsigned int i = 0; i < n_components; ++i)
156 {
157 std::ostringstream namebuf;
158 namebuf << '_' << i;
159 names[i] = name + namebuf.str();
160 }
161 }
162
163 add_data_vector(vec, names);
164}
165
166
167template <int dim, int spacedim>
168template <typename number>
169void
171 const Vector<number> &vec,
172 const std::vector<std::string> &names)
173{
174 Assert(dof_handler != nullptr,
176 // either cell data and one name,
177 // or dof data and n_components names
178 Assert(((vec.size() == dof_handler->get_triangulation().n_active_cells()) &&
179 (names.size() == 1)) ||
180 ((vec.size() == dof_handler->n_dofs()) &&
181 (names.size() == dof_handler->get_fe(0).n_components())),
183 names.size(), dof_handler->get_fe(0).n_components()));
184 for (const auto &name : names)
185 {
186 (void)name;
187 Assert(name.find_first_not_of("abcdefghijklmnopqrstuvwxyz"
188 "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
189 "0123456789_<>()") == std::string::npos,
191 name,
192 name.find_first_not_of("abcdefghijklmnopqrstuvwxyz"
193 "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
194 "0123456789_<>()")));
195 }
196
197 if (vec.size() == dof_handler->n_dofs())
198 {
199 typename std::vector<DataVector>::iterator data_vector = dof_data.begin();
200 for (; data_vector != dof_data.end(); ++data_vector)
201 if (data_vector->names == names)
202 {
203 data_vector->data.reinit(vec.size());
204 std::copy(vec.begin(), vec.end(), data_vector->data.begin());
205 return;
206 }
207
208 // ok. not found. there is a
209 // slight chance that
210 // n_dofs==n_cells, so only
211 // bomb out if the next if
212 // statement will not be run
213 if (dof_handler->n_dofs() !=
214 dof_handler->get_triangulation().n_active_cells())
215 Assert(false, ExcVectorNotDeclared(names[0]));
216 }
217
218 // search cell data
219 if ((vec.size() != dof_handler->n_dofs()) ||
220 (dof_handler->n_dofs() ==
221 dof_handler->get_triangulation().n_active_cells()))
222 {
223 typename std::vector<DataVector>::iterator data_vector =
224 cell_data.begin();
225 for (; data_vector != cell_data.end(); ++data_vector)
226 if (data_vector->names == names)
227 {
228 data_vector->data.reinit(vec.size());
229 std::copy(vec.begin(), vec.end(), data_vector->data.begin());
230 return;
231 }
232 Assert(false, ExcVectorNotDeclared(names[0]));
233 }
234
235 // we have either return or Assert
236 // statements above, so shouldn't
237 // get here!
239}
240
241
242template <int dim, int spacedim>
243void
244DataOutStack<dim, spacedim>::build_patches(const unsigned int nnnn_subdivisions)
245{
246 // this is mostly copied from the
247 // DataOut class
248 unsigned int n_subdivisions =
249 (nnnn_subdivisions != 0) ? nnnn_subdivisions : this->default_subdivisions;
250
251 Assert(n_subdivisions >= 1,
253 n_subdivisions));
254 Assert(dof_handler != nullptr,
256
258
259 const unsigned int n_components = dof_handler->get_fe(0).n_components();
260 const unsigned int n_datasets =
261 dof_data.size() * n_components + cell_data.size();
262
263 // first count the cells we want to
264 // create patches of and make sure
265 // there is enough memory for that
266 unsigned int n_patches = 0;
268 dof_handler->begin_active();
269 cell != dof_handler->end();
270 ++cell)
271 ++n_patches;
272
273
274 // before we start the loop:
275 // create a quadrature rule that
276 // actually has the points on this
277 // patch, and an object that
278 // extracts the data on each
279 // cell to these points
280 QTrapezoid<1> q_trapez;
281 QIterated<dim> patch_points(q_trapez, n_subdivisions);
282
283 // create collection objects from
284 // single quadratures,
285 // and finite elements. if we have
286 // an hp-DoFHandler,
287 // dof_handler.get_fe() returns a
288 // collection of which we do a
289 // shallow copy instead
290 const hp::QCollection<dim> q_collection(patch_points);
291 const hp::FECollection<dim> &fe_collection = dof_handler->get_fe_collection();
292
293 hp::FEValues<dim> x_fe_patch_values(fe_collection,
294 q_collection,
296
297 const unsigned int n_q_points = patch_points.size();
298 std::vector<double> patch_values(n_q_points);
299 std::vector<Vector<double>> patch_values_system(n_q_points,
300 Vector<double>(n_components));
301
302 // add the required number of
303 // patches. first initialize a template
304 // patch with n_q_points (in the plane
305 // of the cells) times n_subdivisions+1 (in
306 // the time direction) points
308 default_patch.n_subdivisions = n_subdivisions;
309 default_patch.reference_cell = ReferenceCells::get_hypercube<dim + 1>();
310 default_patch.data.reinit(n_datasets, n_q_points * (n_subdivisions + 1));
311 patches.insert(patches.end(), n_patches, default_patch);
312
313 // now loop over all cells and
314 // actually create the patches
315 typename std::vector<
317 patches.begin() + (patches.size() - n_patches);
318 unsigned int cell_number = 0;
320 dof_handler->begin_active();
321 cell != dof_handler->end();
322 ++cell, ++patch, ++cell_number)
323 {
324 Assert(cell->is_locally_owned(), ExcNotImplemented());
325
326 Assert(patch != patches.end(), ExcInternalError());
327
328 // first fill in the vertices of the patch
329
330 // Patches are organized such
331 // that the parameter direction
332 // is the last
333 // coordinate. Thus, vertices
334 // are two copies of the space
335 // patch, one at parameter-step
336 // and one at parameter.
337 switch (dim)
338 {
339 case 1:
340 patch->vertices[0] =
341 Point<dim + 1>(cell->vertex(0)[0], parameter - parameter_step);
342 patch->vertices[1] =
343 Point<dim + 1>(cell->vertex(1)[0], parameter - parameter_step);
344 patch->vertices[2] = Point<dim + 1>(cell->vertex(0)[0], parameter);
345 patch->vertices[3] = Point<dim + 1>(cell->vertex(1)[0], parameter);
346 break;
347
348 case 2:
349 patch->vertices[0] = Point<dim + 1>(cell->vertex(0)[0],
350 cell->vertex(0)[1],
352 patch->vertices[1] = Point<dim + 1>(cell->vertex(1)[0],
353 cell->vertex(1)[1],
355 patch->vertices[2] = Point<dim + 1>(cell->vertex(2)[0],
356 cell->vertex(2)[1],
358 patch->vertices[3] = Point<dim + 1>(cell->vertex(3)[0],
359 cell->vertex(3)[1],
361 patch->vertices[4] =
362 Point<dim + 1>(cell->vertex(0)[0], cell->vertex(0)[1], parameter);
363 patch->vertices[5] =
364 Point<dim + 1>(cell->vertex(1)[0], cell->vertex(1)[1], parameter);
365 patch->vertices[6] =
366 Point<dim + 1>(cell->vertex(2)[0], cell->vertex(2)[1], parameter);
367 patch->vertices[7] =
368 Point<dim + 1>(cell->vertex(3)[0], cell->vertex(3)[1], parameter);
369 break;
370
371 default:
373 }
374
375
376 // now fill in the data values.
377 // note that the required order is
378 // with highest coordinate running
379 // fastest, we need to enter each
380 // value (n_subdivisions+1) times
381 // in succession
382 if (n_datasets > 0)
383 {
384 x_fe_patch_values.reinit(cell);
385 const FEValues<dim> &fe_patch_values =
386 x_fe_patch_values.get_present_fe_values();
387
388 // first fill dof_data
389 for (unsigned int dataset = 0; dataset < dof_data.size(); ++dataset)
390 {
391 if (n_components == 1)
392 {
393 fe_patch_values.get_function_values(dof_data[dataset].data,
394 patch_values);
395 for (unsigned int i = 0; i < n_subdivisions + 1; ++i)
396 for (unsigned int q = 0; q < n_q_points; ++q)
397 patch->data(dataset, q + n_q_points * i) =
398 patch_values[q];
399 }
400 else
401 // system of components
402 {
403 fe_patch_values.get_function_values(dof_data[dataset].data,
404 patch_values_system);
405 for (unsigned int component = 0; component < n_components;
406 ++component)
407 for (unsigned int i = 0; i < n_subdivisions + 1; ++i)
408 for (unsigned int q = 0; q < n_q_points; ++q)
409 patch->data(dataset * n_components + component,
410 q + n_q_points * i) =
411 patch_values_system[q](component);
412 }
413 }
414
415 // then do the cell data
416 for (unsigned int dataset = 0; dataset < cell_data.size(); ++dataset)
417 {
418 const double value = cell_data[dataset].data(cell_number);
419 for (unsigned int q = 0; q < n_q_points; ++q)
420 for (unsigned int i = 0; i < n_subdivisions + 1; ++i)
421 patch->data(dataset + dof_data.size() * n_components,
422 q * (n_subdivisions + 1) + i) = value;
423 }
424 }
425 }
426}
427
428
429template <int dim, int spacedim>
430void
432{
433 // release lock on dof handler
434 dof_handler = nullptr;
435 for (typename std::vector<DataVector>::iterator i = dof_data.begin();
436 i != dof_data.end();
437 ++i)
438 i->data.reinit(0);
439
440 for (typename std::vector<DataVector>::iterator i = cell_data.begin();
441 i != cell_data.end();
442 ++i)
443 i->data.reinit(0);
444}
445
446
447
448template <int dim, int spacedim>
449std::size_t
460
461
462
463template <int dim, int spacedim>
464const std::vector<
471
472
473
474template <int dim, int spacedim>
475std::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
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 ReadVector< Number > &fe_function, std::vector< Number > &values) const
Definition point.h:111
unsigned int size() const
virtual size_type size() const override
iterator end()
iterator begin()
const FEValuesType & get_present_fe_values() const
Definition fe_values.h:694
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
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
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.
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
std::enable_if_t< std::is_fundamental_v< T >, 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